matlab解方程，求解惑

[kfzuohh]

在下有一个二阶方程求解，解的方程不难，大概来说就是：
Jr=jr(V1)
Jf=jf(V2)
V1+R*V1/Rsh1+V2+S1*R*Jr-V=0
V2+R*V2/Rsh1+V1+S2*R*Jf-V=0
其中：R=16.7e3;
Rsh1=1e15;
Rsh2=1e15;
r=128e-9;
s1=pi*r*r;
s2=pi*r*r;
V=2.3;
求解V1,V2。
可是我的代码运行的error是：Error: Function definitions are not permitted at the prompt or in scripts.


下面是我的代码：
function y1=jf(V2) %反偏电压V1
%syms V2;
fai2=0.23;%%直接带入数值
q=1.602e-19;
k=1.381e-23;
T=300;
epsil0=9e-12;
pi=3.1416;
h=6.625e-34;
beta=q/(k*T);
m0=9.108e-31;
me=0.023*m0;
A=4*q*pi*me*k^2/h^3;%Richard constant
epsil=10*epsil0;
Nc=2*(2*pi*me*k*T)^1.5/h^3;
Nd=2e17;%%直接带入数值
ek=-log(Nc/Nd)/beta;
E00=13;%%直接带入数值
E0=28;%%直接带入数值
y1=A*T^2*exp(-fai2*beta)*(exp(beta*V2)-1);
end

function y0=jr(V1)%反偏电压V1
%syms V1;
fai1=0.22;%%直接带入数值
q=1.602e-19;
k=1.381e-23;
T=300;
epsil0=9e-12;
pi=3.1416;
h=6.625e-34;
beta=q/(k*T);
m0=9.108e-31;
me=0.023*m0;
A=4*q*pi*me*k^2/h^3;%Richard constant
epsil=10*epsil0;
Nc=2*(2*pi*me*k*T)^1.5/h^3;
Nd=2e17;%%直接带入数值
ek=-log(Nc/Nd)/beta;
E00=13;%%直接带入数值
E0=28;%%直接带入数值
y0=(A*T*(pi*q*E00)^0.5/k*(q*(V1-ek)+q*fai1/((cosh(beta*E00))^2))^0.5*exp(-fai1/E0))*exp(V1*(beta-q/E0));%+A*T^2*exp(-vi*beta)*(1-exp(-V*beta));
end

R=16.7e3;
Rsh1=1e15;
Rsh2=1e15;
r=128e-9;
s1=pi*r*r;
s2=pi*r*r;
V=2.3;
[V1,V2]=solve('V1+R*V1/Rsh1+V2+S1*R*Jr(V1)-V=0','V2+R*V2/Rsh1+V1+S2*R*Jf(V2)-V=0');%需要解的方程组
V1=vap(V1,4);
V2=vap(V2,4);

求大牛给指点一下，到底哪里用错了？？？万分感谢！！