Solving ODEs using trapezoidal method in vector/matrix notation

I have been given an example of solving multiple differential ODEs using the numerical trapezoidal method, and encountering issues when attempting to index vector and matrix values inside a for loop.

Here is the example:

My code is as follows:

clear; clc;

Tend = 0.1; %End time
N = 1000;  %Number of steps 
T = Tend/N; %step size
t = 0:T:Tend; %time span

A = [0 -5.38e6; 1 -47.106]; %System Matrix
B = [5.38e6; 0]; %Input Matrix
I = eye(2);
C1 = inv(I - A.*(T/2));
C2 = (I + A.*(T/2)); 
C3 = C1.*C2; %x(t) coefficient
C4 = (B.*(T/2)); %u(t) coefficient


% Initial Conditions
x1(1) = 0;
x2(1) = 0;

%Create Cell Array
x = zeros(2,1);
xx = size(x);
y = cell(xx);

%Functions
xdot1 = @(x1,t) x1;
xdot2 = @(x2,t) x2;
xdot = @(x1,x2,t)[xdot1;xdot2];

%Assign Values to Cell Array
y{1} = xdot1;
y{2} = xdot2;

%Input
ut = @(t) 187.8.*cos(377.*t);

%Pre-Allocate and store results
xt1 = zeros(2,N+1);

%Run and store values
for n=1:N
    
   xt1(n+1)= C3.*[y{1}; y{2}] + C4.*(ut(t(n)+ut(t(n+1))));
end

I am recieving the following error:

Error using vertcat
Nonscalar arrays of function handles are not allowed; use cell arrays instead.

Any help at finding my mistake would be greatly appreciated

Answer 1

[Chris]

A function handle works like this:

f = @(x,t) x.^2;
y = f(3,12);

the function f takes 3 and 12 as input variables and names them x and t, respectively. Since t is unused in the function, it's irrelevant and may not need to be in the function. It squares x and returns 9. equivalently, f = @(x) x.^2;

A couple things I see:

xdot1 = @(x1,t) x1;

You don't have a function x1 defined--x1 is a vector with an initial value. Calling xdot1(someX,someTime) would always return the entire vector x1 (initially 0). Similarly for xdot2.

xt1(n+1)= C3.*[y{1}; y{2}] + ...

y{1}, at this point, is a function

y{1} = @(x1,t) x1

It can't be concatenated into an array. If you give it numerical input, it will return a number, which is possible to use in an array.

xt1(n+1)= C3.*[y{1}(3,5); y{2})4,2)] + ...

Does this help? I'm sure there will be some more work to do once these errors are cleared up.