# Question:Bug in dsolve numeric or fdiff?

## Question:Bug in dsolve numeric or fdiff?

To summarize,

dsolve numeric parametric form gives wrong answers when fdiff is used to calculate the Hessian. Doing dsolve twice gives the correct answer. Hope I am not making syntax or programming errors.

Also, fdiff is not compatible with vector form.

 > restart;
 > Digits:=15;
 (1)
 > sys:=diff(y1(t),t)=-(u+u^2/2)*y1(t),diff(y2(t),t)=u*y1(t);
 (2)

The system is solved with initial conditions 1,0 and value of u being equal to u1 for t<=0.5, and u2 for t>0.5. The objective is y2 at t =1. There are errors when Hessian is calculated using parametric dsolve and fdiff

 > sol1 := dsolve({sys, y1(0) = alpha, y2(0) = beta}, type = numeric, 'parameters' = [alpha, beta, u],maxfun=0,range=0..0.5):
 > obj1:=proc(u1,u2) local z1,s1,s2,z2; global sol1; sol1('parameters'=[1,0,u1]); z1 := sol1(.5); s1:=subs(z1,y1(t)); s2:=subs(z1,y2(t)); sol1('parameters'=[s1,s2,u2]): z2:=sol1(0.5); -subs(z2,y2(t)); end proc;
 (3)
 > u0:=[0.8,1.8];
 (4)
 > obj1(op(u0));
 (5)
 > Hess1:=Matrix(2,2):for i from 1 to 2 do for j from 1 to 2 do Hess1[i,j]:=fdiff(obj1,[i,j],u0,workprec=1.0);od;od;Hess1;
 (6)

Hessian is wrong. It should be noted that when dsolve is called twice inside the objective function there is no error in Hessian calculation with fdiff. The parametric dsolve is efficient, how to get correct answers with that? Another question is how to use fdiff with vector form instead u1,u2 etc to make it generic for any number of variables?

 > obj2:=proc(u1,u2) local z1,s1,s2,z2,sol2; sol2:=dsolve({op(subs(u=u1,[sys])), y1(0) = 1, y2(0) = 0}, type = numeric,maxfun=0,range=0..0.5); z1 := sol2(.5); s1:=subs(z1,y1(t)); s2:=subs(z1,y2(t)); sol2:=dsolve({op(subs(u=u2,[sys])), y1(0) = s1, y2(0) = s2}, type = numeric,maxfun=0,range=0..0.5): z2:=sol2(0.5); -subs(z2,y2(t)); end proc;
 (7)
 > obj2(op(u0));
 (8)
 > Hess2:=Matrix(2,2):for i from 1 to 2 do for j from 1 to 2 do Hess2[i,j]:=fdiff(obj2,[i,j],u0,workprec=1.0);od;od;Hess2;
 (9)

Doing dsolve twice gives correct answer compared to one parametric dsolve even though objective returns the same number (even gradient returns the expected answers).