Extract of data for the last time step

I have tried to extract theta at the middle of the array for all t. Please, how can I extract theta middle that correspond to the last time step? See the my code:

N = 100;
%% Numerical setup
% step size
h = d/N;

tspan = 0:2000;       % Nonlinear case
nt=length(tspan);

% range of z
z=linspace(0,d,N+1);

% initial conditions
Theta = 0.0001;
theta0 = Theta*sin(pi*z/d);                 % Initial condition: 0
v0 = zeros(1,N+1);
theta0_int=theta0(2:N);
v0_int=v0(2:N);
u0 = [theta0_int'; v0_int'];

% Matrix M
M1=eye(N-1,N-1);
M2=zeros(N-1,N-1);
M=[M1 M2;M2 M2];

% Boundary Conditions
Phi_b = 0;

%% ode solver
counter = 0;
Xis = xi;
for xiidx = 1:length(Xis)
    xi = Xis(xiidx);  
    for Uis = u
        options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
        tic
        [t,y] = ode15s(@(t,y)lcode1(t, y), tspan,u0, options);
        toc
        
        % Extract the solution for theta and v
        theta = [Phi_b*ones(length(t), 1) y(:,1:N-1) Phi_b*ones(length(t), 1)];
        v = [zeros(length(t), 1) y(:,N:(2*N-2)) zeros(length(t), 1)];
        % Extract theta and v data for various xi
        counter = counter + 1;
        THETA(counter,:,:) = theta;
        V(counter,:,:) = v;
        
        % theta at the middle of the layer (i.e., z=d/2)
        theta_middle(counter,:,:) =   theta(:, N/2 +1);
    end
end

[Jan]

I struggle with the term "theta middle at the last for last time step". What does it mean?

[University Glasgow]

I want theta(d/2, maxt(t))

Answer 1

[William Rose]

@University Glasgow, Please post the simplest possible code example that illustrates the issue you have. When you include code, you can include it as actual code, by selecting it and then clicking the "Code" button at the top of the editing window. When you include a function, please also include a (short) script illustrating a call to that function.

You have two nested for loops. Inside the inner loop, you call ode45(). ode45() returns column vector t and array y. Array theta is the same as array y, expcept a column has been added on the right of y and on the left of y. (The added columns have the constant value Phi_b.) You create a 3D array THETA. Each slice of THETA contains the 2D array theta from one pass through the inner loop. THETA is never used, so I'm not sure why it is there. Since tspan is a vector of length 2001 on every pass, t and y will have length 2001 on each pass.

Your code includes the lines

theta_middle(counter,:,:) =   theta(:, N/2 +1);
theta_middle12(counter,:,:) = theta(:, N/2 +1);
theta_middle13(counter,:,:) = theta(:, N/2 +1);
theta_middle14(counter,:,:) = theta(:, N/2 +1);

The lines above save the center column of theta (which is a 2D array) into a slice of theta_middle (which is a 3D array).

Why is theta_middle a 3D array? The right hand side returns a column vector, so you should save these column vectors in a 2D array:

theta_middle(:,counter) =   theta(:, N/2 +1);

The center column of array theta is also the center column of array y. Why not just save the center column of y?

Why do you save the same data the same way in 4 different arrays?

Is it guaranteed that N is even? If N is odd, you will get a non-integer array index in the lines above.

The arrays theta_middle and THETA appear to grow inside a loop. Pre-allocate them for efficiency.

The function Max_Non_linear1_IjuptilK_func() returns vector t. The t vector returned will be from the final loop pass only. This is OK, since t is the same on every pass.

The function Max_Non_linear1_IjuptilK_func() returns z, but no value is assigned to z.

[University Glasgow]

Thank you for this... I actually wanna extract the theta_middle that correspond to the last time step?

[Torsten]

I think you want

% theta at the middle of the layer (i.e., z=d/2)
theta_middle(:,counter) =   theta(:, N/2 +1);

instead of

% theta at the middle of the layer (i.e., z=d/2)
theta_middle(counter,:,:) =   theta(:, N/2 +1);

Then

theta_middle_last_time_step(counter) =  theta_middle(end,counter)

[University Glasgow]

Thank you so much Torsten. This is what i want.