matlab编的龙格库塔法解微分方程，求改错

@[TOmatlab编的龙格库塔法解微分方程，求改错

数学模型：

主程序数据：
R1=0.019; R2=0.022; R3=0.044; R4=0.062;% 一级杆直径，二级杆直径，泵径,油管径
L1=1100; L2=500; %L1碳纤维杆的长度，L2钢杆长度
L=5; %单元节长
h1=1600;h2=1047; %泵挂深度 ，动液面
hc=h1-h2; %泵沉没度
S1=1900; S2=7800; S3=850; %碳纤维的密度，钢杆密度,井液密度
E1=1.110^11;
E2=2.110^11; %钢杆弹性模量
Ps=4.510^6; %沉没压力
Pd=13.510^6; %排出压力
ke=410^6; c0=10;%抽油杆悬挂系统的等效弹簧常数,和阻尼
A1=pi/4R1^2; A2=pi/4R2^2; A3=pi/4R3^2; A4=pi/4*R4^2; %碳纤