どんな数式なのかによってテクニックがあったりなかったりするかと思うんですが。
とりあえず泥臭い方法を書いてみます。
（Matlabにこの辺を扱える関数があったような気がしますが、私は使ったことないです）
ニュートン-ラフソン法を多変数に拡張してやれば解けるんじゃないかと思います。
１変数の時と変わるのは
・微分ではなく偏微分を使う
・割り算が行列の割り算（逆行列を掛ける）になる
ぐらいでしょうか。
４変数は書いてられないので、２変数で。
f1(x1,x2)=0、f2(x1,x2)=0を解くことを考えます。
適当な初期値a,bとして、それと解とのずれをΔx1、Δx2とします。
f1(a+Δx1,b+Δx2)=0、f2(a+Δx1,b+Δx2)=0
となればよいのですが、テイラー展開（だったかな？df1/dx1はf1をx1で微分することを表します）を
使って展開し、高次の微分の項を無視すれば
f1(a+Δx1,b+Δx2)=f1(a,b)+(df1/dx1)(Δx1)+(df1/dx2)(Δx2)
f2(a+Δx1,b+Δx2)=f2(a,b)+(df2/dx1)(Δx1)+(df2/dx2)(Δx2)
となります。
左辺=0として、右辺を行列の形で書くと上式は
[f1(a,b);f2(a,b)]+[df1/dx1,df1/dx2;df2/dx1,df2/dx2]*[Δx1;Δx2]=0
となるので、ここからΔ=[Δx1;Δx2]を求めることができます。
あとはa=a+Δx1、b=b+Δx2としてa,bの初期値を更新し、これを繰り返せば
初期値がいいところにあれば解は収束します。うまくいかなければ発散します。
変数が４つになってもやることは同じです。
Matlabでどう組むかは人によってそれぞれですね。ベタに書くならforループ１つで書けますが、
少しでも楽したいなら最低限f1,f2の関数は作った方がいいでしょう。
私なら関数f1(x1,x2),f2(x1,x2),df11(x1,x2),df12(x1,x2),df21(x1,x2),df22(x1,x2)
を作るかな。あるいはニュートン-ラフソン法以外の方法を探すかも知れません。
他になさそうならあきらめてこれでやるかと。
参考になれば幸いです。
お礼
丁寧な回答ありがとうございました。 ニュートン法自体はしっているのですが、4変数となると・・・「？？？」って感じです（汗） ニュートン法は２次のテイラー展開を行って連立方程式を作り解く方法ですよね？ ニュートン法を4変数に拡張するためには、４変数のテイラー展開を求める必要があるかと思います。 頑張れば導出できるかもしれませんが（自信は全くないですが・・・；ω；）、複雑で長い計算式のためできれば微分をしない方法はないでしょうか？ この方法だと、100～200項くらいは出てきそうなので・・・・。 丁寧なご回答ありがとうございました。
補足
補足： >-（Matlabにこの辺を扱える関数があったような気がしますが、私は使ったことないです） fsolveという関数があり、１変数の非線形方程式の解で、与えられた初期値に最も近い解を導出してくれます。