fit 2nd order differential equation parameters to data

I have a second order differential equation with some tunable parameters that I need to fit to some data. I have spent some time looking over some other postings with answers from (StarStrider especially, with some key assists from Torsten) and have finally gotten my script to run, but something isn't right. I suspect it's something I'm missing about passing in my y(t) data into the function that holds my system of differential equations, but I do not know for sure. I am hoping that someone can see what I'm missing. The data fitted function should follow the data much more closely than it does. The data is y vs t, and the diff. eq. is y''=a, where a is a function of y'. Here is my script:

clear all
close all
clc

Time_data = [ 0.8273 ; 1.1104 ; 1.3940 ; 1.6584 ; 1.8764 ; 2.0173 ; 2.1239 ; 2.2343 ; 2.3498 ]-.35 ; % independent variable, experimental data
Y_data = [ 1.5025 ; 2.4650 ; 3.4274 ; 4.3976 ; 5.3600 ; 6.3302 ; 7.2965 ; 8.2590 ; 9.2253 ] ; % dependent variable, experimental data
a_m_guess = 200 ; %init guess for a_m, the parameter to be fitted
a_m = lsqcurvefit(@(am,t)FittingFunction(am,t),a_m_guess,Time_data,Y_data) ;

figure %plot the data as well as the fitted function
hold 'on'
plot(Time_data, Y_data, 'ko', 'DisplayName','data')
plot_times = linspace(0,max(Time_data),100) ;
plot(plot_times, FittingFunction(a_m,plot_times),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Y'), legend('Location','NorthWest')

function F = FittingFunction(a_m,Time_data)
init_slope = 3.39 ; %This is a boundary condition
tspan = Time_data ;
InitConds = [0 init_slope] ; %y(t=0)=0, y'(t=0)=Us
[t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,a_m),tspan,InitConds) ;
F = y_atTdata(:,1) ; %return y only, not y and y'
end

function Y = SysOfDiffEqns(t,y,a_m)
% y''=a(a_m,y') is the 2nd order diff eq.
% a is a function of a_m() and y'(t)
% need to fit a_m to available pairs of t_data and y_data
% to make system of 1st order deff eqs, 
%    substitute  Y(1) = y     
%                Y(2) = y'
% makes a system of two first order ODEs:
%    Y(1)' = Y(2)
%    Y(2)' = a
Y = zeros(2,1) ; % initialize output matrix
D = 8.7049 ; % constant input parameter
U_m = 0.95*D ; % costant input parameter
U = y(2) ; % U is dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m)) ; %acceleration function, only valid for U<D
a = max(a,0) ; %force to be at least zero
Y(1) = y(2) ; %1st derivative substitution
Y(2) = a ; %2nd derivative equal to a (which is function of y')
end

Thanks very much for any assistance. Again, I'm thinking I somehow need to pass in my y_data into @SysOfDiffEqns and then set y(1)=y_data, but I'm unable to figure it out after several days of fiddling with it.

Answer 1

[Star Strider]

The problem appears to be not allowing the intial conditions to be estimated as parameters as well as the initial value of ‘a_m’.

Also, ‘plot_Times’ must be defined only over the range of times being fitted. If you want to extrapolate to 0 (or any other values), use the deval function.

The only changes I made were to ‘a_m_guess’ in the script, and in the ‘FittingFunction’ function, the ‘InitConds’ assignment and the single parameter ‘a_m(1)’ passed to ‘SysOfDiffEqns’ in the ode45 call.

With only these changes, the code seems to work, and give the appropriate result —

Time_data = [ 0.8273 ; 1.1104 ; 1.3940 ; 1.6584 ; 1.8764 ; 2.0173 ; 2.1239 ; 2.2343 ; 2.3498 ]-.35 ; % independent variable, experimental data
Y_data = [ 1.5025 ; 2.4650 ; 3.4274 ; 4.3976 ; 5.3600 ; 6.3302 ; 7.2965 ; 8.2590 ; 9.2253 ] ; % dependent variable, experimental data
a_m_guess = [200; 1; 1] ; %init guess for a_m, the parameter to be fitted
a_m = lsqcurvefit(@(am,t)FittingFunction(am,t),a_m_guess,Time_data,Y_data) 

Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.

a_m = 3×1

369.4651 1.0532 4.2399

figure %plot the data as well as the fitted function
hold 'on'
plot(Time_data, Y_data, 'ko', 'DisplayName','data')
plot_times = linspace(min(Time_data),max(Time_data),100) ;
plot(plot_times, FittingFunction(a_m,plot_times),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Y'), legend('Location','NorthWest')

function F = FittingFunction(a_m,Time_data)
init_slope = 3.39 ; %This is a boundary condition
tspan = Time_data ;
% InitConds = [0 init_slope] ; %y(t=0)=0, y'(t=0)=Us
InitConds = a_m(2:3);
[t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,a_m(1)),tspan,InitConds) ;
F = y_atTdata(:,1) ; %return y only, not y and y'
end

function Y = SysOfDiffEqns(t,y,a_m)
% y''=a(a_m,y') is the 2nd order diff eq.
% a is a function of a_m() and y'(t)
% need to fit a_m to available pairs of t_data and y_data
% to make system of 1st order deff eqs, 
%    substitute  Y(1) = y     
%                Y(2) = y'
% makes a system of two first order ODEs:
%    Y(1)' = Y(2)
%    Y(2)' = a
Y = zeros(2,1) ; % initialize output matrix
D = 8.7049 ; % constant input parameter
U_m = 0.95*D ; % costant input parameter
U = y(2) ; % U is dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m)) ; %acceleration function, only valid for U<D
a = max(a,0) ; %force to be at least zero
Y(1) = y(2) ; %1st derivative substitution
Y(2) = a ; %2nd derivative equal to a (which is function of y')
end

Make appropriate changes to get different results.

.

[Christopher Neel]

Awesome! Thanks so much. I was looking in the wrong place for the error. It might be a problem to have the initial slope be a fitted parameter- I will see if I can figure out how to force it to a value, but even so this is extremely helpful.

[Star Strider]

As always, my pleasure!

You can use the parameter bounds lsqcurvefit permits as one of its arguments to keep the initial slope relatively tightly constrained, without otherwise changing the code. I would simply constrain it rather than defining it. (The only problem with constraining the parameters is that estimating the confidence intervals of the remaining parameters or of the fit is then not reliable.)

.

[Christopher Neel]

Now it's mostly working, expect that I can't figure out how extrapolate the solution to t=0 as you mentioed in your first comment above. I have two ideas:

1. Use deval, as you suggested. I tried this, as shown in the code below, and it gives me the error "Attempting to evaluate the solution outside the interval [8.273000e-01, 2.349800e+00] where it is defined." If you run the updated script below it will make this error. This is by far what I'd like to do.
2. Do a backwards numerica integration using something like trapz. With the fitted parameters, I have the y'' as a function of y'(t). I think I should be able to start at my first data point of the solved system and numerially integrate backwards. I think this wil be tricky because I'll have to guess the lower limit of y'. Plus it just seems clunky, so I'm hoping not to have to go this route. I haven't tried to code this yet but it's my backup idea.

Here is the improved script, which results in an error from deval. Is there some way to extrapolate the equations to t=0?

clear all
close all
clc

% ydot_o = 3.39 ; %estimate of initial slope
% D_guess = 8.7049 ; %estimate of final slope
Time_data = [ 0.8273 ; 1.1104 ; 1.3940 ; 1.6584 ; 1.8764 ; 2.0173 ; 2.1239 ; 2.2343 ; 2.3498 ] ; % independent variable, experimental data
Y_data = [ 1.5025 ; 2.4650 ; 3.4274 ; 4.3976 ; 5.3600 ; 6.3302 ; 7.2965 ; 8.2590 ; 9.2253 ] ; % dependent variable, experimental data
%note that the data doesn't start at t=0, so the initial conditions are for *the time at which the data starts*
guess = [1000, 1.5, 3.3, 8.7, 0.94]; % [a_m, y(t_o), y'(t_o), D, U_m_factor] init guess for a_m, y_o, and y'_o, D, Um factor
constrain_parameters = 1 ;
if constrain_parameters == 0
    lower_bound = [] ;
    upper_bound = [] ;
elseif constrain_parameters == 1
    lower_bound = [100, 1, 3, 8, 0.91] ;
    upper_bound = [1e4, 2, 4, 9, 0.98] ;
end
return_y_and_ydot = 0 ; % if =0, return only y values (for fitting).  If =1, return both y and y' values (for plotting).
options = optimset() ;
[fitted, resnorm, residuals,~,~] = lsqcurvefit(@(fitted,t)FittingFunction(fitted,t,return_y_and_ydot),guess,Time_data,Y_data,lower_bound,upper_bound,options)

%% evaluation and plots
return_y_and_ydot = 1 ;

figure(1)
subplot(2,1,1) %plot the y(t) data as well as the fitted y(t)
hold 'on'
plot(Time_data, Y_data, 'ko', 'DisplayName','data')
plot_times = linspace(min(Time_data),max(Time_data),1e4) ;
solution = FittingFunction(fitted,plot_times,return_y_and_ydot) ;
plot(plot_times, solution(:,1),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('data and fit vs time'), grid 'on'

subplot(2,1,2) %plot the fitted y'(t) with the transition line
hold 'on'
plot(plot_times, solution(:,2),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Ydot'), legend('Location','NorthWest'), title('ydot vs time'), grid 'on'

figure(2) %plot the residuals
plot(Time_data, residuals,'ko', 'DisplayName','residuals')
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('residuals vs time'), grid 'on'

%% extrapolate solution backwards to t=0
ext_times = linspace(0,max(Time_data),1e4) ; %extend time to zero
ICs = [fitted(2), fitted(3)] ; %what are the new ICs?
sol = ode23s(@(t,y)SysOfDiffEqns(t,y,fitted),Time_data,ICs) ;
ext_y = deval(sol,ext_times) ;
figure(1)
subplot(2,1,1) % the y(t) plot
plot(ext_times,ext_y(1,:),'-g','DisplayName','extension')

%% functions
function F = FittingFunction(guess,Time_data,return_y_and_ydot)
tspan = Time_data ;
InitConds = [guess(2), guess(3)] ; %y(t=t_o), y'(t=t_o)=init_slope
% [t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
[t_atTdata,y_atTdata] = ode23s(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
if return_y_and_ydot == 0 % used during optimization operation against data
    F = y_atTdata(:,1) ; %return y only, not y and y'
elseif return_y_and_ydot == 1 % used to display results
    F = y_atTdata ; %return both y and y'
end
end

function Y = SysOfDiffEqns(t,y,guess)
% y''=a(a_m,y') is the 2nd order diff eq.
% a is a function of a_m() and y'(t)
% need to fit a_m to available pairs of t_data and y_data
% to make system of 1st order deff eqs, 
%    substitute  Y(1) = y     
%                Y(2) = y'
% makes a system of two first order ODEs:
%    Y(1)' = Y(2)
%    Y(2)' = a
Y = zeros(2,1) ; % initialize output matrix
a_m = guess(1) ;
D = guess(4) ; %let parameter be fitted
U_m_factor = guess(5) ; %let parameter be fitted
% D = 8.7049 ; % constant input parameter
% U_m_factor = 0.94 ; % constant input parameter
U_m = U_m_factor*D ; 
U = y(2) ; % U is dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m)) ; %acceleration function, only valid for U<D
a = max(a,0) ; %force to be at least zero
Y(1) = y(2) ; %1st derivative substitution
Y(2) = a ; %2nd derivative equal to a (which is function of y')
end

thanks in advance for any pointers!

[Star Strider]

‘Is there some way to extrapolate the equations to t=0?’

Use the odextend function to calculate values beyond the original region-of-fit.

With respect to deval, my apologies. I cited the wrong function.

.

[Christopher Neel]

That seems to be the perfect function (which somehow I didn't stumble on already today) except that I think it only extends to greater times. I get an error when I try and extend to lower times. I need one that goes the other way. I'm wondering if I have to take my fitted parameters and present it as a new function of ODEs to a boundary value solver within an optimization framework to get what I want. The right hand side boundary would be what I solved for as ICs in my original problem. I would have the optimization solve for the left hand side boundary problems. Is that the best approach?

[Christopher Neel]

I forgot to add what I tried. Here is the text I tried in my last section. It gives the error.

ext_times = linspace(0,max(Time_data),1e4) ; %extend time to zero
ICs = [fitted(2), fitted(3)] ; 
sol = ode23s(@(t,y)SysOfDiffEqns(t,y,fitted),Time_data,ICs) ;
ext_sol = odextend(sol,[],0) ;  %This line gives an error.
ext_y = deval(ext_sol,ext_times) ;

[Christopher Neel]

I figured it out! It finally dawned on me that I can integrate backwards using a regular ode solver by putting in the tspan where it starts at t0 and decreases to zero. That makes all the difference, as it lets me turn my final conditions into initial conditions. In case this helps someone else out in future, here is the full script that does what I wanted. The "extrapolate solution backward" part is the new bit.

clear all
close all
clc

% ydot_o = 3.39 ; %estimate of initial slope
% D_guess = 8.7049 ; %estimate of final slope
Time_data = [ 0.8273 ; 1.1104 ; 1.3940 ; 1.6584 ; 1.8764 ; 2.0173 ; 2.1239 ; 2.2343 ; 2.3498 ] ; % independent variable, experimental data
Y_data = [ 1.5025 ; 2.4650 ; 3.4274 ; 4.3976 ; 5.3600 ; 6.3302 ; 7.2965 ; 8.2590 ; 9.2253 ] ; % dependent variable, experimental data
%note that the data doesn't start at t=0, so the initial conditions are for *the time at which the data starts*
guess = [1000, 1.5, 3.3, 8.7, 0.94]; % [a_m, y(t_o), y'(t_o), D, U_m_factor] init guess for a_m, y_o, and y'_o, D, Um factor
constrain_parameters = 1 ;
if constrain_parameters == 0
    lower_bound = [] ;
    upper_bound = [] ;
elseif constrain_parameters == 1
    lower_bound = [100, 1, 3, 8, 0.91] ;
    upper_bound = [1e4, 2, 4, 9, 0.98] ;
end
return_y_and_ydot = 0 ; % if =0, return only y values (for fitting).  If =1, return both y and y' values (for plotting).
options = optimset() ;
[fitted, resnorm, residuals,~,~] = lsqcurvefit(@(fitted,t)FittingFunction(fitted,t,return_y_and_ydot),guess,Time_data,Y_data,lower_bound,upper_bound,options)

%% evaluation and plots
return_y_and_ydot = 1 ;
trans_threshold = 0.99 ;
D = fitted(4) ;
trans_pt = trans_threshold*D ;

figure(1)

subplot(2,1,1) %plot the y(t) data as well as the fitted y(t)
hold 'on'
plot(Time_data, Y_data, 'ko', 'DisplayName','data')
plot_times = linspace(min(Time_data),max(Time_data),1e4) ;
solution = FittingFunction(fitted,plot_times,return_y_and_ydot) ;
plot(plot_times, solution(:,1),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('data and fit vs time'), grid 'on'

subplot(2,1,2) %plot the fitted y'(t) with the transition line
hold 'on'
plot(plot_times, solution(:,2),'-r','DisplayName','fit') ;
line([min(Time_data),max(Time_data)],[trans_pt, trans_pt],'DisplayName','transition')
xlabel('time'), ylabel('Ydot'), legend('Location','NorthWest'), title('ydot vs time'), grid 'on'

figure(2) %plot the residuals
plot(Time_data, residuals,'ko', 'DisplayName','residuals')
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('residuals vs time'), grid 'on'

%% extrapolate solution backwards to t=0
ext_times = linspace(0,min(Time_data),100) ; %extend time to zero
ICs = [fitted(2), fitted(3)] ; %what are the new ICs?
ext_sol = ode23s(@(t,y)SysOfDiffEqns(t,y,fitted),[min(Time_data), 0],ICs) ;
ext_y = deval(ext_sol,ext_times) ;
figure(1)
subplot(2,1,1) % the y(t) plot
plot(ext_times,ext_y(1,:),'-g','DisplayName','extension')
subplot(2,1,2) % the y'(t) plot
plot(ext_times,ext_y(2,:),'-g','DisplayName','extension')


%% functions
function F = FittingFunction(guess,Time_data,return_y_and_ydot)
tspan = Time_data ;
InitConds = [guess(2), guess(3)] ; %y(t=t_o), y'(t=t_o)=init_slope
% [t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
[t_atTdata,y_atTdata] = ode23s(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
if return_y_and_ydot == 0 % used during optimization operation against data
    F = y_atTdata(:,1) ; %return y only, not y and y'
elseif return_y_and_ydot == 1 % used to display results
    F = y_atTdata ; %return both y and y'
end
end

function Y = SysOfDiffEqns(t,y,guess)
% y''=a(a_m,y') is the 2nd order diff eq.
% a is a function of a_m() and y'(t)
% need to fit a_m to available pairs of t_data and y_data
% to make system of 1st order deff eqs, 
%    substitute  Y(1) = y     
%                Y(2) = y'
% makes a system of two first order ODEs:
%    Y(1)' = Y(2)
%    Y(2)' = a
Y = zeros(2,1) ; % initialize output matrix
a_m = guess(1) ;
D = guess(4) ; %let parameter be fitted
U_m_factor = guess(5) ; %let parameter be fitted
% D = 8.7049 ; % constant input parameter
% U_m_factor = 0.94 ; % constant input parameter
U_m = U_m_factor*D ; 
U = y(2) ; % U is dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m)) ; %acceleration function, only valid for U<D
a = max(a,0) ; %force to be at least zero
Y(1) = y(2) ; %1st derivative substitution
Y(2) = a ; %2nd derivative equal to a (which is function of y')
end

[Star Strider]

I dediced to run the code and post it for the benefit of anyone else who encounters this and wants to see the results —

% ydot_o = 3.39 ; %estimate of initial slope
% D_guess = 8.7049 ; %estimate of final slope
Time_data = [ 0.8273 ; 1.1104 ; 1.3940 ; 1.6584 ; 1.8764 ; 2.0173 ; 2.1239 ; 2.2343 ; 2.3498 ] ; % independent variable, experimental data
Y_data = [ 1.5025 ; 2.4650 ; 3.4274 ; 4.3976 ; 5.3600 ; 6.3302 ; 7.2965 ; 8.2590 ; 9.2253 ] ; % dependent variable, experimental data
%note that the data doesn't start at t=0, so the initial conditions are for *the time at which the data starts*
guess = [1000, 1.5, 3.3, 8.7, 0.94]; % [a_m, y(t_o), y'(t_o), D, U_m_factor] init guess for a_m, y_o, and y'_o, D, Um factor
constrain_parameters = 1 ;
if constrain_parameters == 0
    lower_bound = [] ;
    upper_bound = [] ;
elseif constrain_parameters == 1
    lower_bound = [100, 1, 3, 8, 0.91] ;
    upper_bound = [1e4, 2, 4, 9, 0.98] ;
end
return_y_and_ydot = 0 ; % if =0, return only y values (for fitting).  If =1, return both y and y' values (for plotting).
options = optimset() ;
[fitted, resnorm, residuals,~,~] = lsqcurvefit(@(fitted,t)FittingFunction(fitted,t,return_y_and_ydot),guess,Time_data,Y_data,lower_bound,upper_bound,options)

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

fitted = 1×5

1.0e+03 * 1.1093 0.0015 0.0032 0.0087 0.0009

resnorm = 0.0027

residuals = 9×1

0.0131 -0.0226 -0.0019 0.0205 -0.0099 0.0204 -0.0183 -0.0201 0.0187

%% evaluation and plots
return_y_and_ydot = 1 ;
trans_threshold = 0.99 ;
D = fitted(4) ;
trans_pt = trans_threshold*D ;

figure(1)

subplot(2,1,1) %plot the y(t) data as well as the fitted y(t)
hold 'on'
plot(Time_data, Y_data, 'ko', 'DisplayName','data')
plot_times = linspace(min(Time_data),max(Time_data),1e4) ;
solution = FittingFunction(fitted,plot_times,return_y_and_ydot) ;
plot(plot_times, solution(:,1),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('data and fit vs time'), grid 'on'

subplot(2,1,2) %plot the fitted y'(t) with the transition line
hold 'on'
plot(plot_times, solution(:,2),'-r','DisplayName','fit') ;
line([min(Time_data),max(Time_data)],[trans_pt, trans_pt],'DisplayName','transition')
xlabel('time'), ylabel('Ydot'), legend('Location','NorthWest'), title('ydot vs time'), grid 'on'

figure(2) %plot the residuals
plot(Time_data, residuals,'ko', 'DisplayName','residuals')
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('residuals vs time'), grid 'on'

%% extrapolate solution backwards to t=0
ext_times = linspace(0,min(Time_data),100) ; %extend time to zero
ICs = [fitted(2), fitted(3)] ; %what are the new ICs?
ext_sol = ode23s(@(t,y)SysOfDiffEqns(t,y,fitted),[min(Time_data), 0],ICs) ;
ext_y = deval(ext_sol,ext_times) ;
figure(1)
subplot(2,1,1) % the y(t) plot
plot(ext_times,ext_y(1,:),'-g','DisplayName','extension')
subplot(2,1,2) % the y'(t) plot
plot(ext_times,ext_y(2,:),'-g','DisplayName','extension')

%% functions
function F = FittingFunction(guess,Time_data,return_y_and_ydot)
tspan = Time_data ;
InitConds = [guess(2), guess(3)] ; %y(t=t_o), y'(t=t_o)=init_slope
% [t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
[t_atTdata,y_atTdata] = ode23s(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
if return_y_and_ydot == 0 % used during optimization operation against data
    F = y_atTdata(:,1) ; %return y only, not y and y'
elseif return_y_and_ydot == 1 % used to display results
    F = y_atTdata ; %return both y and y'
end
end

function Y = SysOfDiffEqns(t,y,guess)
% y''=a(a_m,y') is the 2nd order diff eq.
% a is a function of a_m() and y'(t)
% need to fit a_m to available pairs of t_data and y_data
% to make system of 1st order deff eqs, 
%    substitute  Y(1) = y     
%                Y(2) = y'
% makes a system of two first order ODEs:
%    Y(1)' = Y(2)
%    Y(2)' = a
Y = zeros(2,1) ; % initialize output matrix
a_m = guess(1) ;
D = guess(4) ; %let parameter be fitted
U_m_factor = guess(5) ; %let parameter be fitted
% D = 8.7049 ; % constant input parameter
% U_m_factor = 0.94 ; % constant input parameter
U_m = U_m_factor*D ; 
U = y(2) ; % U is dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m)) ; %acceleration function, only valid for U<D
a = max(a,0) ; %force to be at least zero
Y(1) = y(2) ; %1st derivative substitution
Y(2) = a ; %2nd derivative equal to a (which is function of y')
end

I congratulate on extending the original idea!

.