应该可以有toolbox一类的使用，不过我不太常用现成工具。。。
如果自己编程
你可以这样，先把(2)(3)里面的关于fi的项单独提出来变成
比如(2) 定义2(4+8fi)=m, 2(1+fi)^2=n，（方程本身没化成最简形式啊。。）就成了
m(1-2s1)-n/(1-s1)=-1  s1~=1有  m(1-2s1)(1-s1)+1-s1-n=0 展开就可以根据二院一次方程来作出解析解了
(3)目测会出一元三次方程，你可以化成标准形式网上搜下求解公式和判别就行了
然后你看看解出来的是不是有意义，比如虚数解不要
当然可能会得到比较多的解，然后看这架势可能s1 s2都不止一个还要排列组合吧？或者除非你有条件只用某一组特定的解

伪代码大概可以这样

for fi
    call solution_ function2(fi) %% return a vector_s1
    call solution_ function3(fi) %% return a vector_s2
    for 1:size(vector_s1)      for 1:size(vector_s2)
       plot
      end
    end
end

很感兴趣你这个相图画出来会是啥样的。。。

其实比较复杂的就是(3)吧，三元一次，如果实在觉得解析解太麻烦也可以用凑数法，不过不太精确就是了。就是说，从接近正负无穷大开始猜数，用二分法一点点猜过来，算方程的解穿过横轴的位置，程序相对好写，但是比较废机时

或者，你可以直接用Matlab的代数功能来解决这个问题。。
其实我觉得这个方程组的问题，大概可能Maple或者Mathmatic效果会更好点？

多谢了  我还要研究下你提供你的方案，一下没有看明白