Using a for loop for a function inside ODE45

Hello,

I've posted the whole code at the end so that it can make sense if my explantion is a bit wacky.

So I've been trying to find response of spring-mass system with impulse as excitation. The code that i've come up with works, but it requires me write an every inefficient code,

 e = @(time) Y*(exp(lamda.*((t-t1).^2)) + ...
       exp(lamda.*((t-t2).^2)) + ...
       exp(lamda.*((t-t3).^2)) + ...
       exp(lamda.*((t-t4).^2)) + ...
       exp(lamda.*((t-t5).^2)) + ...
       exp(lamda.*((t-t6).^2)) ); 
   edot = @(time) 2*Y*lamda.*((t-t1)*exp(lamda.*((t-t1).^2)) + ...
       (t-t2)*exp(lamda.*((t-t2).^2)) + ...
       (t-t3)*exp(lamda.*((t-t3).^2)) + ...
       (t-t4)*exp(lamda.*((t-t4).^2)) + ...
       (t-t5)*exp(lamda.*((t-t5).^2)) + ...
       (t-t6)*exp(lamda.*((t-t6).^2))); 

I've been trying to find way such that by running a for loop, I don't have to write the code again and again.

Some thing like the one below but I am not able to run this inside a ode function. Please help!

no_cyc = 200;
no_imp = 62;
diff_imp = 0.005;

impulse_timing = (0.005)*[1:1:no_cyc*no_imp];

                %Base excitation
     Y = 0.1; % amplitude
    lamda = -50000000; % exponential fuction 
  


for time_count = tspan
    for impulse_count = [1:1:no_cyc*no_imp]
        e(time_count) = e(time_count) + Y*exp(lambda*(time_count - impulse_timing(imp_count))^2);
        edot(time_count) = edot(time_count) + (2*lambda*(time_count - impulse_timing(imp_count)))*Y*exp(lambda*(time_count - impulse_timing(imp_count)^2));
    end
end

Here is the full code

function twodof_impulse
    % Sampling Parameters
    fs = 2000;  %sampling Frequency
    ts = 1/fs ;     %Sampling interval
    
    %Timespan
    t0 = 0; tf = 40;
    tspan = [t0:ts:tf];
    
    
    %initial conditions
    y1_0 = 0; v1_0 =0 ;
    y2_0 = 0; v2_0 = 0;
    
    IC = [y1_0,v1_0, y2_0, v2_0];
    
    
    %Numerical Integration and Assignment of variables
    [time, state_values] = ode45(@g,tspan,IC);
    y1 = state_values(:,1);
    v2 = state_values(:,2);
    y2 = state_values(:,3);
    v2 = state_values(:,4);
   
    nfft = length(y1); %length of Time domain signal
    nfft2 = 2^nextpow2(nfft); %length of the signal in power of 2
    grid_points = 2*nfft2;
    
    %Base excitation
     Y = 0.1; % amplitude
    lamda = -100; % exponential fuction 
   t1 = 5;   % Time of impulses
    t2 = 10;
    t3 = 15;
    t4 = 20;
    t5 = 25;
    t6 = 30;

    
   e = Y*(exp(lamda.*((tspan-t1).^2)) + ...
       exp(lamda.*((tspan-t2).^2)) + ...
       exp(lamda.*((tspan-t3).^2)) + ...
       exp(lamda.*((tspan-t4).^2)) + ...
       exp(lamda.*((tspan-t5).^2)) + ...
       exp(lamda.*((tspan-t6).^2)) ); 
    
    
  %FFT for y1 and y2
    
    ff1 = fft(y1, grid_points);
    fff1 = ff1(1:nfft2/2);
    
    ff2 = fft(y2, grid_points);
    fff2 = ff2(1:nfft2/2);
    
    f = fs*(0:nfft2/2 - 1)/grid_points; %Conversion to from parts to Hz
    
    %Graph Controls
    
    xlength_linear = nfft2*((1*3*3)/(10*2*5*2*3));
    xlength_log = nfft2*((1*3*3)/(10*2*5*2*3));
    %ploting results
    
    %ploting Input
    figure(1)
    clf
    
     subplot(311)
    plot(tspan,e), xlabel('Time(s)'),ylabel('Displacement(m)')
    title('Base Excitaion')
    
    %first plot of y1(Wheel)
    subplot(312)
    plot(time,y1), xlabel('Time(s)'),ylabel('Displacement(m)')
    title('Wheel Response')
    
    %second plot of y2(Body)
    subplot(313)
    plot(time,y2), xlabel('Time(s)'), ylabel('Displacement(m)')
    title('Body Response')
    
    
figure(2)

    %ploting frequncy analysis y1(Wheel)
    subplot(411)
    plot(f(1:xlength_linear), abs(fff1(1:xlength_linear)))
    title ("Frequency vs amplitude of wheel displacement (Linear)");
    xlabel('Frequency,Hz');
    ylabel('Amplitude,m^2');

    
    %ploting frequncy analysis y2(Body)
    subplot(413)
    plot(f(1:xlength_linear), abs(fff2(1:xlength_linear)))
    
    title ("Frequency vs Amplitude of body dispacement (Linear)");
    xlabel('Frequency,Hz');
    ylabel('Amplitude,m^2');

%ploting y1(semilog)
    subplot(412)
    semilogy(f(1:xlength_log) , abs(fff1(1:xlength_log)))
    
    title ("Frequency vs Amplitude of wheel dispacement (Log)");
    xlabel('Frequency,Hz');
    ylabel('Amplitude,m^2');

    %ploting y2(semilog)
    subplot(414)
    semilogy(f(1:xlength_log) , abs(fff2(1:xlength_log)))
    
    title ("Frequency vs Amplitude of body dispacement (Log)");
    xlabel('Frequency,Hz');
    ylabel('Amplitude,m^2');


end

%Ode 45 fuction
function sdot = g(t,s)
    
    %system Varibles(All in SI Units)
    Mw = 35;Mb = 200;  
    k1 = 180000; k2 = 17000;
    c1 = 10; c2 = 200;
   
    %force in the base( Here Wheel )
    Y = 0.1; % amplitude
    lamda = -100; % exponential fuction 
    t1 = 5;   % Time of impulses
    t2 = 10;
    t3 = 15;
    t4 = 20;
    t5 = 25;
    t6 = 30;

    
   e = @(time) Y*(exp(lamda.*((t-t1).^2)) + ...
       exp(lamda.*((t-t2).^2)) + ...
       exp(lamda.*((t-t3).^2)) + ...
       exp(lamda.*((t-t4).^2)) + ...
       exp(lamda.*((t-t5).^2)) + ...
       exp(lamda.*((t-t6).^2)) ); 
   edot = @(time) 2*Y*lamda.*((t-t1)*exp(lamda.*((t-t1).^2)) + ...
       (t-t2)*exp(lamda.*((t-t2).^2)) + ...
       (t-t3)*exp(lamda.*((t-t3).^2)) + ...
       (t-t4)*exp(lamda.*((t-t4).^2)) + ...
       (t-t5)*exp(lamda.*((t-t5).^2)) + ...
       (t-t6)*exp(lamda.*((t-t6).^2))); 
   
   %State function
   
   sdot = [s(2);
       (-1/Mw)*((c1+c2)*s(2) + -c2*s(4) + (k1 + k2)*s(1) - k2*s(3) -(c1*edot(t) + k1*e(t)));
       s(4);
       (-1/Mb)*((-c2)*s(2) + c2*s(4) + (-k2)*s(1) + k2*s(3)) ];
   
   
   end

[Stephen23]

Both of the anonymous functions e and edot have one input argunment time, but this input is completely ignored.

"I've been trying to find way such that by running a for loop, I don't have to write the code again and again."

At first glance it looks like you are trying to do the ODE solver's job for it. The whole point of ODE solvers is that they estimate the ODE behavior by repeatedly calling the input function with suitable input values. It is normally not required to loop over the tspan as you are trying to do, that is what the ODE solver does.

Why do you think that you need those loops inside the function g ?

[Ram Hembram]

e = @(time) Y*(exp(lamda.*((t-t1).^2)) + ...
       exp(lamda.*((t-t2).^2)) + ...
       exp(lamda.*((t-t3).^2)) + ...
       exp(lamda.*((t-t4).^2)) + ...
       exp(lamda.*((t-t5).^2)) + ...
       exp(lamda.*((t-t6).^2)) );
edot = @(time) 2*Y*lamda.*((t-t1)*exp(lamda.*((t-t1).^2)) + ...
       (t-t2)*exp(lamda.*((t-t2).^2)) + ...
       (t-t3)*exp(lamda.*((t-t3).^2)) + ...
       (t-t4)*exp(lamda.*((t-t4).^2)) + ...
       (t-t5)*exp(lamda.*((t-t5).^2)) + ...
       (t-t6)*exp(lamda.*((t-t6).^2)));

"Why do you think that you need those loops inside the function g ?"

Here I'm trying to give Impulse as an input. So at t1,...,t6 the function gives kind of a spike and for all other values it goes to almost zero. Here since its only six its no problem to write, but suppose I want give 1000 impulses then this code becomes reduntant and long. So inside the g function(ODE45) I was trying find way by using a for loop or anything else ( here to make the funtion e and e dot).

I was trying to find a better way to do this. Below is what I've been doing If I wanted 20 impulses.

e = @(time) Y*(exp(lamda.*((t - impulse_timing(1)).^2)) + ... 
exp(lamda.*((t-impulse_timing(2)).^2)) + ...
exp(lamda.*((t-impulse_timing(3)).^2)) + ...
exp(lamda.*((t-impulse_timing(4)).^2)) + ...
exp(lamda.*((t-impulse_timing(5)).^2)) + ...
exp(lamda.*((t-impulse_timing(6)).^2)) + ...
exp(lamda.*((t-impulse_timing(7)).^2)) + ...
exp(lamda.*((t-impulse_timing(8)).^2)) + ...
exp(lamda.*((t-impulse_timing(9)).^2)) + ...
exp(lamda.*((t-impulse_timing(10)).^2)) + ...
exp(lamda.*((t-impulse_timing(11)).^2)) + ...
exp(lamda.*((t-impulse_timing(12)).^2)) + ...
exp(lamda.*((t-impulse_timing(13)).^2)) + ...
exp(lamda.*((t-impulse_timing(14)).^2)) + ...
exp(lamda.*((t-impulse_timing(15)).^2)) + ...
exp(lamda.*((t-impulse_timing(16)).^2)) + ...
exp(lamda.*((t-impulse_timing(17)).^2)) + ...
exp(lamda.*((t-impulse_timing(18)).^2)) + ...
exp(lamda.*((t-impulse_timing(19)).^2)) + ...
exp(lamda.*((t-impulse_timing(20)).^2))) 

[Stephen23]

Copy-and-pasting code like that is a sign that you are doing something wrong. The simple and efficient MATLAB approach would be to use operations on the entire vectors/arrays, all wrapped in sum, e.g. assuming that t and lamda are scalar:

Y*sum(exp(lamda.*(t-impulse_timing).^2))

As far as I can tell, there is no need to define this in an anonymous function.

Why do you define an input argument (to the anonymous function) which is not used?