You are over-writing the symboic variables, psi and theta by taking input from the user. Rename the variables which you will use for substitution
Also, the if-else condition block seems redundant, as regardless of the value of phi, psi = -phi
syms phi theta psi Psi
h = 110.75;
d2r = pi/180;
r2d = 180/pi;
ROT = [(cos(phi)*cos(psi))-(cos(theta)*sin(phi)*sin(psi)) -(cos(phi)*sin(psi))-(cos(theta)*cos(psi)*sin(phi)) (sin(phi)*sin(theta));...
(cos(psi)*sin(phi))+(cos(phi)*cos(theta)*sin(psi)) (cos(phi)*cos(theta)*cos(psi)-sin(phi)*sin(psi)) -(cos(phi)*sin(theta));...
(sin(theta)*sin(psi)) (cos(psi)*sin(theta)) cos(theta)];
px = (h/2)*(ROT(1,1)-ROT(2,2));
py = -h*ROT(2,1);
eqn = ROT(1,2) == ROT(2,1);
Psi = eqn;
pz = 435; %input('Please enter value for pz: ');
phi0 = 5; %input('Please enter value for phi: ')*d2r;
theta0 = 10; %input('Please enter value for theta: ')*d2r;
if phi0 >= 0
psi0 = -phi0 ;
elseif phi0 <= 0
psi0 = -phi0 ;
end
jr = [diff(px, pz), diff(px, phi), diff(px, theta); diff(py, pz), diff(py, phi), diff(py, theta); 1, 0, 0; 0, 1, 0; 0, 0, 1; diff(Psi, pz), diff(Psi, phi), diff(Psi, theta)];
Jr = double(subs(jr,{phi,theta,psi},{phi0,theta0,psi0}))
