The MATLAB and theoretical solutions are equivalent.
Matlab is giving you this:
d/dz( cot(z) ) = ( cot(z)^2 + 1 )
Let's work from there...
= ( 1/tan(z)^2 + 1 ) % By definition cot(z)=1/tan(z)
= ( cos(z)^2 / sin(z)^2 + 1 ) % Using tan(z) = sin(z)/cos(z)
= ( (1  sin(z)^2) / sin(z)^2 + 1 ) % Using cos(z)^2 + sin(z)^2 = 1
= ( 1/sin(z)^2  sin(z)^2/sin(z)^2 + 1 ) % Expanding the fraction
= ( 1/sin(z)^2  1 + 1 ) % Simplifying
= 1/sin(z)^2 % Expected result!
So the MATLAB result is exactly as expected from theory, as you stated at the end. d is treated as a constant in the partial derivative, so remains as a coefficient.
If you didn't fancy the manual workingout, you could just get MATLAB to check the equivalence for you...
simplify( R_z  (d/sin(z)^2) ) % = 0, so they are the same.
You can get a the expected form from MATLAB, by using rewrite and simplify .
% rewrite the result R_z in terms in the SIN function
simplify(rewrite( R_z, 'sin' ) )
% >> ans = d/sin(z)^2
However, note that this is not well defined. From the simplify docs
Simplification of mathematical expression is not a clearly defined subject. There is no universal idea as to which form of an expression is simplest.
You are better off using the simplify(...) = 0 approach shown above, and then printing the desired result to screen if you want to. If you're not outputting the result string then the expression's form is irrelevant and you can march on regardless!
As for wolframalpha, I cannot reproduce your issue...
But again, you can note that the first result you show is equivalent:
% Your result from wolfram alpha
syms x
diff(cot(x), x)
% =  cot(x)^2  1
% =  ( cot(x)^2 + 1 )
% This is the same as the MATLAB result, same reasoning follows
