Help to solve a large system of Nonlinear equations

I tried to solve a system of 41 nonlinear equation - 41 variables.
Physically,This system describe a very hyperstatically concrete girder
behavior under large loads.

I have try to solve it with a software "Systems of Nonlinear
Equations" found on www.numericalmathematics.com. This software use
generalized Newton method, but this software can't take a system of 41
equations.

I have try to solve it with mathematics library PETSc but this
algortihm need that I build the Jacobian matrix of 41 by 41 element.
It's very hard because I have about 400 element different of zero.
Probably, some software build Jacobian matrix.

Finally, the difficulty of this system is to found initial value of
the 41 independant variables for that it converge.

The system is:

41 variables:
P,M0,M1,M2,M3,M4,Vba,Vca,Vda,Vea,f1p,f2p,f3p,f4p,K1,K3,K4,K5,K6,K7,K8,Lam1,Lam2,
Lam3,Lam4,f1c,f2c,f3c,f4c,c1,c2,c3,c4,T1x,T2x,T3x,T4x,Tde,Tcd,Tbc,Tab

41 equations:
"-M0+L*P/4-L/8*tan(f4c)*T4x-L/8*tan(f3c)*(T3x+T4x)-L/8*tan(f2c)*(T2x+T3x+T4x)-L/8*tan(f1c)*(T1x+T2x+T3x+T4x)+T4x*Vea+T3x*(Vea-Vba)+T2x*(Vea-Vca)+T1x*(Vea-Vda)-T4x*e-T3x*e-T2x*e-T1x*e=0"
"-M1+3*L*P/16-L/8*tan(f4c)*T4x-L/8*tan(f3c)*(T3x+T4x)-L/8*tan(f2c)*(T2x+T3x+T4x)+T4x*Vda+T3x*(Vda-Vba)+T2x*(Vda-Vca)-T4x*e-T3x*e-T2x*e-T1x*e/2=0"
"-M2+L*P/8-L/8*tan(f4c)*T4x-L/8*tan(f3c)*(T3x+T4x)+T4x*Vca+T3x*(Vca-Vba)-T4x*e-T3x*e-T2x*e/2=0"
"-M3+L*P/16-L/8*tan(f4c)*T4x+T4x*Vba-T4x*e-T3x*e/2=0"
"-M4-T4x*e=0"
"Vba-K1*sin(Lam4*L/8)+e*cos(Lam4*L/8)-L/8*tan(f4c)-e+P*L/(16*T4x)=0"
"Vca-K3*sin(Lam3*L/4)-K4*cos(Lam3*L/4)-L/8*tan(f3c)-e-1/(T3x+T4x)*(-P*L/8+L/8*tan(f4c)*T4x+T3x*Vba)=0"
"Vda-K5*sin(3*Lam2*L/8)-K6*cos(3*Lam2*L/8)-L/8*tan(f2c)-e-1/(T2x+T3x+T4x)*(-3*P*L/16+L/8*tan(f4c)*T4x+L/8*tan(f3c)*T3x+L/8*tan(f3c)*T4x+T3x*Vba+T2x*Vca)=0"
"Vea-K7*sin(Lam1*L/2)-K8*cos(Lam1*L/2)-L/8*tan(f1c)-e-1/(T1x+T2x+T3x+T4x)*(-P*L/4+L/8*tan(f4c)*T4x+L/8*tan(f3c)*(T3x+T4x)+L/8*tan(f2c)*(T2x+T3x+T4x)+T3x*Vba+T2x*Vca+T1x*Vda)=0"
"f1p-K5*Lam2*cos(3*Lam2*L/8)+K6*Lam2*sin(3*Lam2*L/8)-tan(f2c)+P/(2*(T2x+T3x+T4x))=0"
"f2p-K3*Lam3*cos(Lam3*L/4)+K4*Lam3*sin(Lam3*L/4)-tan(f3c)+P/(2*(T3x+T4x))=0"
"f3p-K1*Lam4*cos(Lam4*L/8)-e*Lam4*sin(Lam4*L/8)-tan(f4c)+P/(2*T4x)=0"
"f4p-K1*Lam4-tan(f4c)+P/(2*T4x)=0"
"K1*sin(Lam4*L/8)-e*cos(Lam4*L/8)-(T3x+T4x)/T4x*(K3*sin(Lam3*L/8)+K4*cos(Lam3*L/8))-T3x/T4x*e=0"
"K1*Lam4*cos(Lam4*L/8)+e*Lam4*sin(Lam4*L/8)+tan(f4c)-P/(2*T4x)-K3*Lam3*cos(Lam3*L/8)+K4*Lam3*sin(Lam3*L/8)-tan(f3c)+P/(2*(T3x+T4x))=0"
"K3*sin(Lam3*L/4)+K4*cos(Lam3*L/4)-(T2x+T3x+T4x)/(T3x+T4x)*(K5*sin(Lam2*L/4)+K6*cos(Lam2*L/4))-T2x/(T3x+T4x)*e=0"
"K3*Lam3*cos(Lam3*L/4)-K4*Lam3*sin(Lam3*L/4)+tan(f3c)-P/(2*(T3x+T4x))-K5*Lam2*cos(Lam2*L/4)+K6*Lam2*sin(Lam2*L/4)-tan(f2c)+P/(2*(T2x+T3x+T4x))=0"
"K5*sin(3*Lam2*L/8)+K6*cos(3*Lam2*L/8)-(T1x+T2x+T3x+T4x)/(T2x+T3x+T4x)*(K7*sin(3*Lam1*L/8)+K8*cos(3*Lam1*L/8))-T1x/(T2x+T3x+T4x)*e=0"
"K5*Lam2*cos(3*Lam2*L/8)-K6*Lam2*sin(3*Lam2*L/8)+tan(f2c)-P/(2*(T2x+T3x+T4x))-K7*Lam1*cos(3*Lam1*L/8)+K8*Lam1*sin(3*Lam1*L/8)-tan(f1c)+P/(2*(T1x+T2x+T3x+T4x))=0"
"K7*Lam1*cos(Lam1*L/2)-K8*Lam1*sin(Lam1*L/2)-P/(2*(T1x+T2x+T3x+T4x))+tan(f1c)=0"
"Lam1-sqr((T1x+T2x+T3x+T4x)/(Ec*I))=0"
"Lam2-sqr((T2x+T3x+T4x)/(Ec*I))=0"
"Lam3-sqr((T3x+T4x)/(Ec*I))=0"
"Lam4-sqr(T4x/(Ec*I))=0"
"tan(f1c)-(Vea-Vda+(dt-c1)*(1-cos(f1p)))/(L/8+(dt-c1)*sin(f1p))=0"
"tan(f2c)-(Vda-Vca-(dt-c1)*(1-cos(f1p))+(dt-c2)*(1-cos(f2p)))/(L/8-(dt-c1)*sin(f1p)+(dt-c2)*sin(f2p))=0"
"tan(f3c)-(Vca-Vba-(dt-c2)*(1-cos(f2p))+(dt-c3)*(1-cos(f3p)))/(L/8-(dt-c2)*sin(f2p)+(dt-c3)*sin(f3p))=0"
"tan(f4c)-(Vba-(dt-c3)*(1-cos(f3p))+(dt-c4)*(1-cos(f4p)))/(L/8-(dt-c3)*sin(f3p)+(dt-c4)*sin(f4p))=0"
"c1-Yg-(T1x+T2x+T3x+T4x)*Ec*I/(M1*(Aac*Eac+Ac*Ec))=0"
"c2-Yg-(T2x+T3x+T4x)*Ec*I/(M2*(Aac*Eac+Ac*Ec))=0"
"c3-Yg-(T3x+T4x)*Ec*I/(M3*(Aac*Eac+Ac*Ec))=0"
"c4-Yg-T4x*Ec*I/(M4*(Aac*Eac+Ac*Ec))=0"
"T1x-cos(f1c)*Et*At1*Tde=0"
"T2x-cos(f2c)*Et*At2*(Tde+Tcd)/2=0"
"T3x-cos(f3c)*Et*At3*(Tde+Tcd+Tbc)/3=0"
"T4x-cos(f4c)*Et*At4*(Tde+Tcd+Tbc+Tab)/4=0"
"Tde-(L/8+(dt-c1)*sin(f1p))/(L/8*cos(f1c))+1=0"
"Tcd-(L/8-(dt-c1)*sin(f1p)+(dt-c2)*sin(f2p))/(L/8*cos(f2c))+1=0"
"Tbc-(L/8-(dt-c2)*sin(f2p)+(dt-c3)*sin(f3p))/(L/8*cos(f3c))+1=0"
"Tab-(L/8-(dt-c3)*sin(f3p)+(dt-c4)*sin(f4p))/(L/8*cos(f4c))+1=0"
"Def0+M0*Yg/(I*Ec)+1/(Aac*Eac+Ac*Ec)*(T1x+T2x+T3x+T4x)=0"
 
15 Constants: dt=650 L=20000 Yg=288.53 e=361.47 I=5546495774 Ac=105000
              Aac=10000 At1=2000 At2=2000 At3=2000 At4=2000 Ec=43866
              Eac=200000 Et=200000 Def0=-0.0030


I need help about the method to use, a convenient software, and a
method to find initial value of the variables so the system converge ?

Sincerely