how to importing m-file with internal function to simulink?

I have a program of simulating a "train model" with a main file by the name 'main.m' with some codes. in 'main.m' file I used ode45 to solve differential equations of Newton Law of train model.

as you know, with ode45 we need a function to describe the differential equations and for this, I used another file by the name of 'vdp.m'.

Now I wanna use 'main.m' file program as a block in Simulink, because this codes are describing my model behavior.

Is there any way to import 'main.m' file to a block and use the block within my simulink model or a way to combining both two m-files into one m-file and making a function for using MATLAB Function Block?

below there is my files.

'main.m':

[t,x]=ode45(@vdp,[0 100],[0 0 -15 0]);
subplot(2,1,1);
plot(t,x(:,1),'r-',t,x(:,3),'b-');
title('Positions');
legend('Loco','Wagon');
xlabel('Time');
ylabel('Distance');
grid;

'vdp.m':

function dx = vdp(t,x)
%%Setting Parameters
c0 = 7.6658*10^-3; % unit = Nkg^-1
cv = 1.08*10^-4;   % unit = Ns(mkg)^-1
ca = 2.06*10^-5;   % unit = Ns^2(m^2 kg)^-1
m1 = 50000;       % unit = kg
m2 = 48500;       % unit = kg
k1 = 85*10^2;      % unit = Nm^-1
d1 = 85*10^4;      % unit = kgs^-1
% Force Input
u = 3000;       % unit = N  48750
%u2 = 0;
teta1 = 0;
teta2 = 0;
D1 = 0;
D2 = 0;
% dx=zeros(6,1); % a column vector
%%Generate Control Input
if t>=0 && t<=10
    a=0;
    u1=u*a;
elseif t>10 && t<=15
    a=1;
    u1=u*a;
elseif t>15 && t<=55
    a=50;
    u1=u*a;
elseif t>55 && t<=75
    a=-97;
    u1=u*a;
else
    a=0;
    u1=u*a;
end
%%State Equations
% x1 -> position of loco
% x2 -> velocity of loco
% x3 -> position of wagon
% x4 -> velocity of wagon

if t>=0 && t<=10
    dx(1)=0; dx(2)=0; dx(3)=0; dx(4)=0;
    dx = [dx(1);dx(2);dx(3);dx(4)];
elseif t>10 && t<=75
    dx(1)=x(2);
    dx(2)=(1/m1)*(u1-k1*(x(1)-x(3))-d1*(x(2)-x(4))-(c0+cv*x(2))*m1- ...
        ca*((x(2))^2)*(m1+m2)-9.98*(sin(teta1))*m1-0.004*D1*m1);
    dx(3)=x(4);
    dx(4)=(1/m2)*(-k1*(x(3)-x(1))-d1*(x(4)-x(2))-(c0+cv*x(4))*m2- ... 
        9.98*(sin(teta2))*m2-0.004*D2*m2);
    dx = [dx(1);dx(2);dx(3);dx(4)];
else
    dx(1)=0; dx(2)=0; dx(3)=0; dx(4)=0;
    dx = [dx(1);dx(2);dx(3);dx(4)];
end