Question: Differental equations (problem with dsolve ?)

Hi,

I have to simulate some kinematic movements of machine. There are tree differential equations (ce2,ce2 amd ce3) which describe movement (there are angular velocity ac1, ac2, ac3 or angular beta, ang, teta ). I know exact solution in MatLab but I am trying do it the same in Maple (13). Unfortunatelly results-curves I recieved are not the same (In Matlab I have very smooth). Is there any Maple expert who can check my equations and give me some suggestions to solve my problem, please? I do not know what I am doing wrongly.  (Please, find exact solution-curves below).
I would like to mention that initial conditions (IC) for example beta has a huge impact on curves. For IC=0 I am not receiving curves (only for teta). If I change IC for beta (for example to 0.0000001 or 0.01) I am receiving someting but it is not I wolud like. But anyway, I should recieve good results for intial condition = 0. Where is the problem?
 
Many thanks for your help!

> restart;
> with(linalg, subvector, submatrix, inverse, extend, copyinto, diag, transpose);
> with(LinearAlgebra);
> with(plots);
>
> a1 := l[0]*sin(beta(t))+l[3]*sin(beta(t)+ang(t));
              
> a2 := l[3]*sin(beta(t)+ang(t));
                       
> a3 := R[0]*cos(teta(t));
                             

> b1 := -l[0]*cos(beta(t))-l[3]*cos(beta(t)+ang(t));
              
> b2 := -l[3]*cos(beta(t)+ang(t));
> b3 := R[0]*sin(teta(t));
                             
> c1 := 1;
                                      1
> c2 := 0;
                                      0
> c3 := 0;
                                      0
> J := Array([[a1, a2, a3], [b1, b2, b3], [c1, c2, c3]]);
                           Array(%id = 152599672)
> MatrixInverse(`<|>`(`<,>`(a1, a2, a3), `<,>`(b1, b2, b3), `<,>`(c1, c2, c3)));
  
#ja2, ja3, jb2, jb3, jc2, jc3 - there are values from inverse matrix                        

> ja2 := sin(teta(t))/(l[3]*(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t))));
          (sin(teta(t)))/(l[3] (cos(teta(t)) cos(beta(t) + ang(t)) + sin(beta(t) + ang(t)) sin(teta(t))))
> ja3 := cos(beta(t)+ang(t))/(R[0]*(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t))));
      (cos(beta(t) + ang(t)))/(R[0] (cos(teta(t)) cos(beta(t) + ang(t) + sin(beta(t) + ang(t)) sin(teta(t))))
> jb2 := -cos(teta(t))/(l[3]*(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t))));
         - (cos(teta(t)))/(l[3] (cos(teta(t)) cos(beta(t) + ang(t)) + sin(beta(t) + ang(t)) sin(teta(t))))
> jb3 := sin(beta(t)+ang(t))/(R[0]*(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t))));
      (sin(beta(t) + ang(t)))/(R[0] (cos(teta(t)) cos(beta(t) + ang(t)) + sin(beta(t) + ang(t)) sin(teta(t))))
> jc2 := -(sin(teta(t))*l[0]*sin(beta(t))+sin(teta(t))*l[3]*sin(beta(t)+ang(t))+cos(teta(t))*l[0]*cos(beta(t))+cos(teta(t))*l[3]*cos(beta(t)+ang(t)))/(l[3]*(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t))));
- (sin(teta(t)) l[0] sin(beta(t)) + sin(teta(t)) l[3] sin(beta(t) + ang(t)) + cos(teta(t)) l[0] cos(beta(t)) + cos(teta(t)) l[3] cos(beta(t) + ang(t)))/ (l[3] (cos(teta(t)) cos(beta(t) + ang(t)) + sin(beta(t) + ang(t)) sin(teta(t))))
> jc3 := -l[0]*(cos(beta(t)+ang(t))*sin(beta(t))-sin(beta(t)+ang(t))*cos(beta(t)))/(R[0]*(cos(teta(t))*cos(beta(t)+ang(t))+sin(beta(t)+ang(t))*sin(teta(t))));
  - (l[0] (cos(beta(t) + ang(t)) sin(beta(t)) - sin(beta(t) + ang(t)) cos(beta(t))))/(R[0] (cos(teta(t)) cos(beta(t) + ang(t)) + sin(beta(t) + ang(t)) sin(teta(t))))

> ce1 := -(diff(beta(t), t))+.5*sin(beta(t))*ja2+.5*tan(beta(t))*ja3 = 0;
> ce2 := -(diff(ang(t), t))+.5*sin(beta(t))*jb2+.5*tan(beta(t))*jb3 = 0;
> ce3 := -(diff(teta(t), t))+.5*sin(beta(t))*jc2+.5*tan(beta(t))*jc3+.5*cos(beta(t)) = 0;

> F := solve({ce1, ce2, ce3}, {diff(ang(t), t), diff(beta(t), t), diff(teta(t), t)});

> ac1 := unapply(subs(F, diff(beta(t), t)), t);
> ac2 := unapply(subs(F, diff(ang(t), t)), t);
> ac3 := unapply(subs(F, diff(teta(t), t)), t);
 

 
> l[0] := 1; l[1] := .5; l[3] := .35; R[1] := .15; R[0] := 2; tw := 10; gama := 1.5; T := 20; i[4] := 10.2;
                                   
                              
> ini := beta(0) = 0.1e-1, ang(0) = 0., teta(0) = 0.;
                 
> G := dsolve({ce1, ce2, ce3, ini}, {ang(t), beta(t), teta(t)}, type = numeric, method = classical[foreuler], output = listprocedure, stepsize = 0.1e-2);
> cb := odeplot(G, [t, ac1(t)], 0 .. 10, color = black);
> cbb := odeplot(G, [t, ac2(t)], 0 .. 10, color = red);
> cbbb := odeplot(G, [t, ac3(t)], 0 .. 10, color = blue);
> display([cb], [cbb], [cbb]);

Curves shold looked like this: