How to simplify this and make it efficient?

lc, clear all

A= 0.06;
k_l = 26400; %Linear stiffness
m = 483; %Mass
l =0.5; 
d =-0.005;
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
Om_array = linspace(0,20,20); %in rad/s-1
l_array = linspace(0,1,20);
[om_array, L_array] = meshgrid(Om_array, l_array);
Response_amp = zeros(size(Om_array));
T = 150;
x0 = [0,0];
for i=1:numel(Om_array)
    for j=1:numel(l_array)
        Om = om_array(i,j);
        l = L_array(i,j);
        k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
        
        f = @(t,x) [ x(2); ...
        -(2*k_s*(x(1)-(A*sin(Om*t)))).* ...
        (sqrt((l-d).^2 + (x(1)-(A*sin(Om*t)))^2) - l)/ ...
        (m*(sqrt((l-d).^2 + (x(1)-(A*sin(Om*t)))^2))) ];
        
        [t, x] = ode45(f,[100,T],x0);
        Response_amp(i,j) = (max(x(:,1)) - min(x(:,1)))/2;
        %    xval(i) = Om/(2*pi) ;
    end
end
%% plot
figure(1);
ax = axes();
view(3);
hold(ax);
view([30 33]);
grid on
mesh(om_array/(2*pi),L_array,Response_amp) ;
xlabel('Frequency (Hz)')
ylabel('Length of the spring (m)')
zlabel('Response Amplitude (m)')
set(gca,'FontSize',15)
% set(gca,'xtick',[])
% set(gca,'ytick',[])
% set(gca,'ztick',[])

%%
%l = linspace(0,1,40);
%b = max(max(Response_amp));

hold on
d =-0.01;
Om_array = linspace(0,20,20); %in rad/s-1
l_array = linspace(0,1,20);
[om_array, L_array] = meshgrid(Om_array, l_array);
Response_amp = zeros(size(Om_array));
T = 150;
x0 = [0,0];
for i=1:numel(Om_array)
    for j=1:numel(l_array)
        Om = om_array(i,j);
        l = L_array(i,j);
        k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
        
        f = @(t,x) [ x(2); ...
        -(2*k_s*(x(1)-(A*sin(Om*t)))).* ...
        (sqrt((l-d).^2 + (x(1)-(A*sin(Om*t)))^2) - l)/ ...
        (m*(sqrt((l-d).^2 + (x(1)-(A*sin(Om*t)))^2))) ];
        
        [t, x] = ode45(f,[100,T],x0);
        Response_amp(i,j) = (max(x(:,1)) - min(x(:,1)))/2;
        %    xval(i) = Om/(2*pi) ;
    end
end
%% plot
figure(1);
ax = axes();
view(3);
hold(ax);
view([30 33]);
grid on
mesh(om_array/(2*pi),L_array,Response_amp) ;
xlabel('Frequency (Hz)')
ylabel('Length of the spring (m)')
zlabel('Response Amplitude (m)')
set(gca,'FontSize',15)

mesh(om_array/(2*pi),L_array,Response_amp) ;

hold off

Hi, all. This code shows 2 different graphs in the same figure when d = -0.005 and -0.01. However, I wish to simplify this code as it seems to be quite long-winded. Also, I want to vary d from -0.005 to -0.03.

Thanks for reading and I would appreciate your time for solving this problem.