Dynamic size differential equation system

I would like to create a system of N differential equations in a function and solve them with ode45. I tried to create them in a for loop:

function output = deqns(par)
alpha=par.alpha;
beta=par.beta;
n=par.n; % n = N/2 where N is the number of equations

ydot=zeros(2*n,1);
for k = 1:2*n
    if k <= n
        ydot(k) = alpha*(optV(par,y(k+n))-y(k))+beta*(y(mod(k+n,n)+1)-y(k));
    else
        ydot(k) = y(mod(k+n,n)+1)-y(k);   
    end      
end
output = ydot;
end

The result should look like this, where the total number of equations is 2*n. (This is my original system)

When I try to solve this with ode45, i get a bunch of errors, because 'y' is not defined in the For loop. I tried to solve the problem but all I achived was some different errors.

[t, y] = ode45(@(t,y) deqns(par), [t0,tmax], init);

ERRORS:
Undefined function 'y' for input arguments of type 'double'.
Error in deqns (line 10)
        ydot(k) = alpha*(optV(par,y(k+n))-y(k))+beta*(y(mod(k+n,n)+1)-y(k));
Error in asd>@(t,y)deqns(par)
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 asd (line 49)
[t, y] = ode45(@(t,y) deqns(par), [t0,tmax], init);

It is important that I need to avoid symbolic calculations (if possible). Any help is appreciated.

Answer 1

[Star Strider]

I have absolutely no idea what you’re doing. Regardless, your ‘deqns’ function must have as two of its arguments your independent variable (here apprently ‘t’) and your dependent variable (here apparently ‘y’).

With those changes, you function runs without error, however you will probably want to change it so it runs much more efficiently.

Try this:

function output = deqns(t,y,par)
alpha=par.alpha;
beta=par.beta;
n=par.n; % n = N/2 where N is the number of equations
ydot=zeros(2*n,1);
for k = 1:2*n
    if k <= n
        ydot(k) = alpha*(optV(par,y(k+n))-y(k))+beta*(y(mod(k+n,n)+1)-y(k));
    else
        ydot(k) = y(mod(k+n,n)+1)-y(k);   
    end      
end
output = ydot;
end


optV = @(b1,b2) 42;                                             % Insert Correct Function
par.alpha = rand;                                               % Assign Correct Value
par.beta = rand;                                                % Assign Correct Value
par.n = 4;                                                      % Assign Correct Value

t0 = 0;                                                         % Assign Correct Value
tmax = 1;                                                       % Assign Correct Value
init = zeros(2*par.n,1);                                        % Assign Correct Value
[t, y] = ode45(@(t,y) deqns(t,y,par), [t0,tmax], init);

figure
plot(t, y)
grid

Experiment to get the result you want.

[Zsolt Fischer]

Thank you, adding t and y as arguments solved the problem!

[Star Strider]

As always, my pleasure!

Answer 2

[Steven Lord]

There's no need for a loop here. Work on vectors.

function dvhdt = myodefun(t, vh, alpha, optV, beta)
% Compute N from the vh vector with which the ODE solver calls your function
N = numel(vh)/2;

% Unpack v and h into separate vectors
v = vh(1:N);
h = vh(N+1:end);

% Compute the updated dv/dt and dh/dt separately
% Use the computation for dh/dt to update dv/dt
dhdt = [diff(v); v(1)-v(N)];
dvdt = alpha*(optV - v)+beta*dhdt;

% Pack dv/dt and dh/dt together
dvhdt = [dvdt; dhdt];
end

You will need to pass alpha, optV, and beta into your ODE function as extra parameters.