Errors using ODE45: "Output of the function must be the same size as the input"

I'm trying to solve a 2nd order ODE numerically where the ODE includes a double integral that can't be solved analytically. The integral is a function of 3 variables: y, z1 and z2. I would like to integrate over z1 and z2 to get a function of y (I used 'int' not 'integral' so I can leave the y in place, but put in integration bounds for z1 and z2 so it can go back and evaluate for a particular value of y), then use ODE45 to solve my equation for y.

My code is as follows:

g3 = @(y) int(int((exp(-z1.^2).*exp(-z2.^2)/(2.*pi)).*tanh(sqrt(c0 - y.^2/c0).*z1)...
    .* tanh(sqrt(c0).*z2),z2, -inf, inf), z1, -inf, inf);

[V] = odeToVectorField(diff(y, 2) == y - h.^2 .* g3(y) - 2.*n.*kroneckerDelta(t));

M = matlabFunction(V,'vars', {'t','Y'});

[T, Y] = ode45(M,[0 1],[c0 0]);

So g3 is just expressing my integral as a function of the variable (y) that is eventually the subject of my ODE. z1 and z2 are dummy integration variables. c0, h, n are all scalars.

The error I get is as follows:

Error using integralCalc/finalInputChecks (line 515)
Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the
'ArrayValued' option to true.

Error in integralCalc/iterateScalarValued (line 315)
                finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 103)
        [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in
symengine>@(t,Y)[Y(2);-integral(@(z1)integral(@(z2)(exp(-z1.^2).*exp(-z2.^2).*tanh(z1.*sqrt(Y(1).^2.*(-2.749426546736005e1)+3.637122079828446e-2)).*tanh(z2.*1.907124033677004e-1).*(1.0./2.0))./pi,-Inf,Inf),-Inf,Inf)+Y(1)-(t==0.0).*(1.0./5.0)]

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in random_pot_v2 (line 44)
[T, Y] = ode45(M,[0 1],[c0 0]);

And I have no idea whether the error is because I've expressed the integral wrongly and/or because I've done something wrong with the ODE45. I know that solving the ODE numerically is possible because I've seen the result elsewhere. Any help much appreciated.

Answer 1

[Torsten]

Try

M=@(t,Y)[Y(2); -integral(@(z1)integral(@(z2)(exp(-z1.^2).*exp(-z2.^2).*tanh(z1.*sqrt(Y(1).^2.*(-2.749426546736005e1)+3.637122079828446e-2)).*tanh(z2.*1.907124033677004e-1).*(1.0./2.0))./pi,-Inf,Inf,'ArrayValued',true),-Inf,Inf,'ArrayValued',true)+Y(1)-(t==0.0).*(1.0./5.0)];
[T, Y] = ode45(M,[0 1],[c0 0]);

And the "kroneckerDelta(t)" is superfluous ; it won't influence the solution (at least numerically), I guess.

Best wishes

Torsten.

[CM]

Thank you for taking the time to reply. Agreed re the delta function.

I've run the code and although it's now error-free, it's giving me a trajectory that diverges at large times. Inside the integral, am I right in thinking that the Y(1) inside the square root is now just a scalar constant (equal to my initial condition), or does it act as part of the ODE that is solved in the ODE45 expression?

[Torsten]

Y(1) inside the function expression for M is the time-dependent solution of the ODE at time t.

By the way: Your nested integral is separable.

You should test whether

M=@(t,Y)[Y(2);-0.5*pi*integral(@(z1)(exp(-z1.^2).*tanh(z1*sqrt(Y(1)^2*(-2.749426546736005e1)+3.637122079828446e-2))),-Inf,Inf)*integral(@(z2)(exp(-z2.^2).*tanh(z2*1.907124033677004e-1)),-Inf,Inf)+Y(1)-0.2*(t==0.0)];

works faster.

You could even evaluate the second integral in advance since it doesn't depend on Y.

Best wishes

Torsten.