求助！！！哪位大侠能帮帮忙看一下！求，高人！！！

天一阁

最初的方程组如下：
f1(θ3,θ4)=r2cosθ2+r3cosθ3-r1-r4cosθ4=0
f2(θ3, θ4)=r2sinθ2+r3sinθ3-r4sinθ4=0求助！！！！求助！！！哪位大侠能给出以下程序的牛顿——辛普森方程求解的M文件。
最终能绘制出一下程序的图形来！！！非常感谢！！！！！！！！
% 设定各杆的长度(单位:毫米)
rs(1)=304.8;      
rs(2)=101.6;        
rs(3)=254.0;   
rs(4)=177.8;        
dr=pi/180.0;
th(1)=0.0;         
th(2)=45*dr;      
th(3)=135*dr;      
th(3)=pi-asin(sin(th(2))*rs(3)/rs(4))
dth=5*dr;                              
% 曲柄输入角从0度变化到360度,步长为5度,计算th34
for i=1:72
  [th3,th4]=ntrps(th,rs);
    th34(i,:)=[th(1)/dr th3/dr th4/dr];
    th(1)=th(1)+dth;                  
    th(2)=th3;                        
    th(3)=th4;                          
end
% 绘制输出角th(2)与th(3)—输入角th(1)的关系曲线
subplot(2,2,1)                        
plot(th34(:,1),th34(:,2),th34(:,1),th34(:,3))
axis([0 360 0 170])
grid                                 
ylabel('从动件角位移/deg')
title('角位移线图')
text(110,110,'摇杆4角位移')
text(50,35,'连杆3角位移')