How to get fmincon to work with an ode45?

Hello fellow members of the community

I've written a code for the purpose of finding the specific initial values of an ODE of thrown particle. I wanted it be based on ode45 and fmincon function.

So much to say, it is not working and spits out an error

I'll be very thankful for any help

Cheers

Jan

PS

Here I've copied full code, that should work on any computer.

fun = @Throw; 
x0 = [10 6];
nonlcon = @Ceilling;
[v0opt, fval] = fmincon(fun,x0,[],[],[],[],[],[],nonlcon);

function dpos = Throw(t,pos)  
mu = 0.0094;
g = 10;
dpos = zeros(4,1);
dpos(1) = pos(3);
dpos(2) = pos(4);
dpos(3) = -mu * pos(3)* sqrt(pos(3)^2 + pos(4)^2);
dpos(4) = -g - mu * pos(4) * sqrt(pos(3)^2 + pos(4)^2);
end

function j = FindParam(v0)
targetPosZ = 5;
targetPosY = 2;
Opt = odeset('Events', @StopIntegration);
[t,pos] = ode45(@Throw,[0, 6],[0, 1, v0(1), v0(2)],Opt);
plane = plot(pos(:,1),pos(:,2),[targetPosZ, targetPosZ],[0 3.5],targetPosZ,targetPosY,'o');
j = (targetPosY - pos(end,2))^2 % obj. fun. = (dist. to target)^2
end

function [c,ceq] = Ceilling(v0)
ceq = [];
targetPosZ = 5;
targetPosY = 2;
Opt = odeset('Events', @StopIntegration);
[t,pos] = ode45(@Throw,[0, 10],[0, 1, v0(1), v0(2)],Opt);
h = v0(2)^2/2*g;
v = sqrt(v0(2)^2 + v0(1)^2);
c(1) = max(h - 5);    % height of the ceilling
c(2) = min(y(2) + 0); % projectile can't touch the floor
c(3) = max(v - 20);   % max velocity of 20
end

function [value, isterminal, direction] = StopIntegration(t, y)
targetPosZ = 5;              % integration stops when projectile hits the wall 
value = (y(1) > targetPosZ); 
isterminal = 1;
direction  = 0;
end

Answer 1

[Star Strider]

See: Optimizing a Simulation or Ordinary Differential Equation

[JV]

I found it - I had used wrong function handle. Still, I will for sure look at that documentation, so thanks!

[Star Strider]

As always, my pleasure!