Solving system of ODEs in Matlab

Question

I am trying to code an ODE system for $n=[-2,-1,0,1]$

$\frac{dx_n}{dt}=V_nx_n-x_{n+1}-x_{n-1}$

with some IC, say $[0 ~ 0~ 0~0]$: I have also a B.C at the two ends $n=-2$ and $n=1$, which I don't know how to implement at the moment, so to start I just tried to code with only I.Cs but having some problem. The code is

tic;
x0 = [0 0 0 0];
tspan = [0, 2];
V=-2.5;

eoms = @(t,x) [V*x(1)-x(2); V*x(2)-x(3)-x(1); V*x(3)-x(4)-x(2); V*x(4)-x(5)-x(3)];
[t, x] = ode45(eoms, tspan, x0)
toc
plot(t,x)

where x(1) corresponds to $x_{-2}$ , x(2) corresponds to $x_{-1}$ and so on, for the four equations corresponding to $n=-2,-1,0,1$. The term $x_{-3}$ has been ignored. The code gives some errors "Index exceeds matrix dimensions."

Answer 1

function testODE
k=pi/2; p=[-40:40];
TT=1;
z=zeros(1,numel(p)-2);
w=-2*cos(k);
gamma=1;
z=zeros(1,numel(p)-2)
v=-2.5;
v1=-2.625;
v2=-2.375;
psi3=TT*exp(3*i*k);
psi2=TT*exp(2*i*k);
psi1=-psi3+(v2-w+gamma*abs(psi2).^2).*psi2;
psi0=-psi2+(v1-w+gamma*abs(psi1).^2).*psi1;
R0=(psi0.*exp(-i*k)-psi1)./(exp(-i*k)-exp(i*k));
R1=(psi0.*exp(i*k)-psi1)./(exp(i*k)-exp(-i*k));
y = (p <= 1) .* (R0 * exp(1i*k*p) + R1 * exp(-1i*k*p)) + (p >= 2) .* (TT * exp(1i*k*p));
z(1+(numel(p)-2)/2:2+(numel(p)-2)/2)=[v*y(1+ceil(end/2)),v*y(2+ceil(end/2))]
f = @(t,y) [(R0 * exp(1i*k*p(1)) + R1 * exp(-1i*k*p(1))).*exp(-1i*w*t); diff(y,2)+2*y(2:end-1); (TT * exp(1i*k*p(end))).*exp(-1i*w*t)];
y0 = y;
t0 = 0;
T = 5;
N = 50;
[Y,t] = RK42d(f, t0, T, y0, N);
figure;
hold on;
plot(p,Y,'r');
plot(p,Y,'b');
end

function [Y, t] = RK42d(f, t0, T, y0, N)
%the input of f and y0 must both be vectors
%composed of two fuctions and their respective
%intial conditions
%the vector Y will return a vector as well
h = (T - t0)/(N - 1); %Calculate and store the step size
Y = zeros(N,81); %Initialize the X and Y vector
t = linspace(t0,T,N); % A vector to store the time values
Y(1,:) = y0'; % Start Y vector at the intial values.
for i = 1:(N-1)
    y = Y(i,:)';
    k1 = f(t(i),y);
    k2 = f(t(i) +0.5*h, y +0.5*h*k1);
    k3 = f(t(i) +0.5*h, y +0.5*h*k2);
    k4 = f(t(i) +h    , y +h*k3);
    Y(i+1,:) = y + (h/6)*(k1+ 2.*k2 + 2*k3 + k4);
    %Update approximation
end
end