how to save only the final output of ode 45 solve

%% RSFODE-CRACK_III QUASI-STATIC APPROXIMATION
% Parameter values as in Rubin & Ampuero (2005)
tic
V_init=10e-9;                                   %in meters/s
Tau_dot_inf=10e-2;                         % Pa/s  remote stressing rate 
sigma=100;                                      %Mpa normal stress 
mu_dash=11.56;                             %Gpa rigidity modulus. 
options = odeset('RelTol', 1e-35, 'AbsTol', 1e-6, 'MaxStep', 1e5);
Gmod = 3e4;                                % in MPa
v_shear = 3e9;                             % in microns/sec
eta= 0.5*Gmod/v_shear;                     % Radiation damping term

%RATE AND STATE PARAMETERS

D_c=1;                                     %in mm
b=0.01;                                    % rate state parameter that characterize the magnitude of evolution effect.
a_b=0.5;                                   % a/b ratio
a=a_b*b;                                   % Rate state parameter that characterize the magnitude of direct effect.


%FAULT LENGTH AND DISCRETIZATION ELEMENT
Nx=2^12;                                   % No. of discretization element
spacing=0.05;                              % spacing between the points 
Lx=Nx*spacing;                             % Length of the fault in x-direction.

% Creating a random distribution of v values 
V_init_dist=10e-12 + (10e-11-10e-12).*rand(Nx,1);
%kernel 
Kx = 2*pi/Nx * [0:Nx/2-1, 0, -Nx/2+1:-1];
%parameters that need to be passed on to the derivate script
par = [a,b,D_c,sigma,Tau_dot_inf,mu_dash];
law = 'A';
%create time intervals for the iteration
time_span_interval=linspace(0,1e9,2e7);
t_initial=0;                               % Initial time 
% Create a distribution of theta 
theta_init=(D_c./V_init_dist);


% PREALLOCATING a single vector for initial conditions 
initial_conditions=zeros(2*Nx,1) ;
initial_conditions(1:Nx,1)=theta_init;
initial_conditions(Nx+1:2*Nx,1)=V_init_dist;

%SOLVING FOR V AND THETA 

for i=1:3%50%length(time_span_interval)
    t_i=time_span_interval(i);
    t_f=time_span_interval(i+1);
    
    [t,y] = ode45(@(t,y)rsfode_crack_III(t,y,Nx,par,law),...
        [0 (t_f-t_i)/2 (t_f-t_i)],initial_conditions,options);
    
    initial_conditions=y(end,:);
    
    
end 

Unrecognized function or variable 'rsfode_crack_III'.

Error in solution>@(t,y)rsfode_crack_III(t,y,Nx,par,law) (line 51)
[t,y] = ode45(@(t,y)rsfode_crack_III(t,y,Nx,par,law),...

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

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Answer 1

[Sam Chak]

Hi @SOUVIK DARIPA

I created a dummy model to test the for-loop code.

Though you didn't annotate the code properly, I believe that time span should be

[t_i, t_f]

instead of

[0 (t_f-t_i)/2 (t_f-t_i)]

%% RSFODE-CRACK_III QUASI-STATIC APPROXIMATION
% Parameter values as in Rubin & Ampuero (2005)
tic
V_init=10e-9;                              %in meters/s
Tau_dot_inf=10e-2;                         % Pa/s  remote stressing rate 
sigma=100;                                 %Mpa normal stress 
mu_dash=11.56;                             %Gpa rigidity modulus. 
options = odeset('RelTol', 3e-14, 'AbsTol', 1e-6, 'MaxStep', 1e5);
Gmod = 3e4;                                % in MPa
v_shear = 3e9;                             % in microns/sec
eta= 0.5*Gmod/v_shear;                     % Radiation damping term

%RATE AND STATE PARAMETERS

D_c=1;                                     %in mm
b=0.01;                                    % rate state parameter that characterize the magnitude of evolution effect.
a_b=0.5;                                   % a/b ratio
a=a_b*b;                                   % Rate state parameter that characterize the magnitude of direct effect.

%FAULT LENGTH AND DISCRETIZATION ELEMENT
% Nx=2^12;                                   % No. of discretization element
Nx=3;
spacing=0.05;                              % spacing between the points 
Lx=Nx*spacing;                             % Length of the fault in x-direction.

% Creating a random distribution of v values 
V_init_dist=10e-12 + (10e-11-10e-12).*rand(Nx,1);
%kernel 
Kx = 2*pi/Nx * [0:Nx/2-1, 0, -Nx/2+1:-1];
%parameters that need to be passed on to the derivate script
par = [a,b,D_c,sigma,Tau_dot_inf,mu_dash];
law = 'A';

% create time intervals for the iteration
% time_span_interval=linspace(0,1e9,2e7);
time_span_interval=linspace(0, 60, 61);
t_initial=0;                               % Initial time 

% Create a distribution of theta 
theta_init=(D_c./V_init_dist);

% PREALLOCATING a single vector for initial conditions 
initial_conditions              = zeros(2*Nx, 1);
initial_conditions(1:Nx,1)      = theta_init;
initial_conditions(Nx+1:2*Nx,1) = V_init_dist;

% SOLVING FOR V AND THETA 
for i = 1:length(time_span_interval)-1
    t_i   = time_span_interval(i);
    t_f   = time_span_interval(i+1);
    [t,y] = ode45(@(t,y) rsfode_crack_III(t, y, Nx, par, law), [t_i t_f], initial_conditions, options);
    initial_conditions=y(end,:);
    plot(t, y), hold on
end 

% [t,y] = ode45(@(t,y) rsfode_crack_III(t,y,Nx,par,law), time_span_interval, initial_conditions, options);
% plot(t, y), hold on
grid on

%% Dummy ODE model
function dy = rsfode_crack_III(t,y,Nx,par,law)
    A   = [zeros(2*Nx-1, 1), eye(2*Nx-1);
           - par];
    B   = [zeros(2*Nx-1, 1); 1];
    K   = lqr(A, B, eye(2*Nx), 1);
    u   = - K*y;
    dy  = A*y + B*u;
end