it seems that your boundary condition y(0) causes second derivative to be indefinite (~NaN) diff(y(x),x)|(x=0) = -(x/y(x))(x=0) = -(0/y(0))=-(0/0)=indefinite. Also y(x) in the denominator causes singularities along the solution range for a few sample initial conditions i tried. You can also try numerical solutions via following fragment > sol := dsolve({de,y(1)=1,D(y)(1)=0},numeric,range=0..2); > sol(0.1); > plots[odeplot](sol);

You can convert those trigonometric equations to polynomial form by applying substitution cos(theta)=(1-t^2)/(1+t^2) and sin(t)=2*t/(1+t^2) (with t=tan(theta/2)) for each theta_i, and after this conversion (to polynomial form) apply Grobner basis to solve this new polynomial system. But this process can be also not too fast or not memory efficient. Note: In kinematic analysis of industrial manipulator, people deal with this complexity by spherical-wrist orientation/position decoupling if present/applicable.