自己编译的求解微分方程m文件不知能否供VB输入参数调用DLL

[jiangnanpg]

自己编译的m文件不知能否供VB输入参数调用DLL进行运算并在VB中输出结果，其中L、W、H、p、D为VB的5个文本框输入参数，THETA_A为VB的一个文本框输出参数
（m文件运算结果），求大神指导（在MATLAB中给这5个参数赋值，可以得出结果）
function [THETA_A,zdlyl]=bxyjs(L,W,H,p,D)
clear;close;clc;format long
global THETA
global THETA_A
global L
global W
global H
global p
global D
THETA=0.8;
THETA_A=2;
while abs(THETA_A-THETA)>0.001
THETA=THETA+0.0001;
y0=0;
[x1,y1]=ode45(@f,[0,L],y0);
THETA_A=max(y1);
end
t1=9.160178682401550e+008;
[x2,y2]=ode45(@q,[0:0.001:THETA_A],t1);
zdlyl=max(y2);
plot(x2,y2)
grid on
xlabel('变形角θ/rad'),ylabel('应力σmax/MPa')
title('不同压强下的应力变化曲线')
function dydx=f(x,y)% dθ/ds
global THETA
global W
global H
global p
global D
E=20.5e10;%梁的弹性模量
a=0.021;
b=0.033;%铰链，座头尺寸
dydx=(2*(pi*D^2*p/4)*(cos(THETA/2)-cos(y-THETA/2))/(E*(W*H^3/12))+((pi*D^2*p/4)*(b*cos(THETA/2)+a*sin(THETA/2))/(E*(W*H^3/12)))^2)^0.5;
function dndm=q(m,n)% dσmax/dθ
global THETA_A
global W
global H
global p
global D
E=20.5e10;%梁的弹性模量
a=0.021;
b=0.033;%铰链，座头尺寸
dndm=6*(pi*D^2*p/4)*sin(m-THETA_A/2)/(W*H^2)/((2*(pi*D^2*p/4)*(cos(THETA_A/2)-cos(m-THETA_A/2))/(E*(W*H^3/12))+((pi*D^2*p/4)*(b*cos(THETA_A/2)+a*sin(THETA_A/2))/(E*(W*H^3/12)))^2)^0.5)+(pi*D^2*p/4)*sin(THETA_A/2-m)/(W*H);