关于Matlab求解含变量的二次方程组，程序有什么问题？

Question

需要对四个方程进行求解后(解出来n1,n3是关于J的解)，再进行画图(以J为横坐标，E为纵坐标画图)，需求如附图所示。自己写的程序如下：
syms n1 n2
[n1,n2]=solve('(3.2E7+3.1E7)n1==5.6E3n2-2E-11n1n2+0.255E-15n2^2+J/4*(1.5E-31.6E-19)','3.1E7n1==(5.6E3+8.2E2)n2+1.255E-15n2^2-3J/(41.5E-31.6E-19)','n1','n2');
syms n3 n4
[n3,n4]=solve('(3.2E7+3.1E7)n3=5.6E3n4+J/(41.5E-31.6E-19)','3.1E7n3=(5.6E3+8.2E2)n4-3J/(41.5E-3*1.6E-19)','n3','n4');

J=(0:100);
E=13.6*n1/n3;
plot(J,E)

麻烦各位帮忙指导指导，万分谢谢！
文献里面是仿真出下图的红实线，实在搞不懂，来论坛请教各位彦祖了！

Answer 1

solve函数得到的是解析表达式，然后可以用subs()函数给参数赋值，所以solve表达式里的参数不需要赋具体值，这样不容易出错。第二组方程n3,n4只有1组解，第一组方程n1,n2有3组解，这里我只选了第一组也就是n1(1)和n3(1)。

syms n1 n2;
[n1,n2]=solve('(k1+k2)*n1==k3*n2-k4*n1*n2+a*k5*n2^2+J/(4*d*e)','k2*n1==(k3+k6)*n2+(1+a)*k5*n2^2-3*J/(4*d*e)','n1','n2');
syms n3 n4;
[n3,n4]=solve('(k1+k2)*n3=k3*n4+J/(4*d*e)','k2*n3=(k3+k6)*n4-3*J/(4*d*e)','n3','n4');

k1 = 3.2E7;
k2 = 3.1E7;
k3 = 5.6E3;
k4 = 2E-11;
k5 = 5E-15;
k6 = 8.2E2;
a  = 0.25;
d  = 1.5E-3;
e  = 1.6E-19;

%J=0.01:1:100;
J = logspace(-2,2,20);
E=16.3*n1(1)/n3(1);

loglog(J,subs(E),'-s');
ylim([0.8 25]);
yticks([1 5 10 20]);
ylabel('\etaEQE (%)');
xlim([0.01 100]);
xticklabels({'0.01','0.1','1','10','100'});
xlabel('Current Density (mA/cm^2)');
set(gca, 'box','off');