help improving my code or debugging it for Galerkin method

Question

I'm trying to plot it but I cannot, where it keeps running without giving any plots for more than 40 mins, and I'm quite sure about code...

The issue that I don't know to how to debug the code at all, so please debug it for me if u can or fix my code if it has an issue for u

(*Geometrical data*)
d = 2*10^-6;  h = 2*10^-6;  b = 2*10^-6; l = 200*10^-6; EE = 1.66*10^11;  Eps = 8.854*10^-12;  ρ = 2332;

(*Numerical problem data*)
Subscript[α, 2] = (6*Eps*l^4)/(EE*h^3*
d^3); Subscript[V, pull] = Sqrt[1.72/Subscript[α, 2]]; Subscript[V, dc] = (Range[100] - 1)*Subscript[V, pull]/100;
Solut = {};

(*Number of modes*)
n = 3; 

(*Normalized frequencies*)
Subscript[ω, 1] = 3.51602; Subscript[ω, 2] = 22.0345; \Subscript[ω, 3] = 61.6972; Subscript[σ, 1] = 0.7341; \Subscript[σ, 2] = 1.0185; Subscript[σ, 3] = 0.9992;

(*Modeshapes*)
Subscript[ϕ, 1][x_] = Cosh[Sqrt[Subscript[ω, 1]] x] - Cos[Sqrt[Subscript[ω, 1]] x] - Subscript[σ, 
1]*(Sinh[Sqrt[Subscript[ω, 1]] x] - 
  Sin[Sqrt[Subscript[ω, 1]] x]);
Subscript[ϕ, 2][x_] = Cosh[Sqrt[Subscript[ω, 2]] x] - Cos[Sqrt[Subscript[ω, 2]] x] - Subscript[σ, 
2]*(Sinh[Sqrt[Subscript[ω, 2]] x] - 
  Sin[Sqrt[Subscript[ω, 2]] x]);
Subscript[ϕ, 3][x_] = Cosh[Sqrt[Subscript[ω, 3]] x] - Cos[Sqrt[Subscript[ω, 3]] x] - Subscript[σ, 
3]*(Sinh[Sqrt[Subscript[ω, 3]] x] - 
  Sin[Sqrt[Subscript[ω, 3]] x]);

(*The derived eq for Cantilver beam *)
SysEq1 = (\!\(\*SubsuperscriptBox[\(∫\), \(0\), \(1\)]\(\(\*SubscriptBox[\(ϕ\), \(j\)]\)[x] \*SuperscriptBox[\((1 - \*UnderoverscriptBox[\(∑\), \(l = 1\), \(n\)]\*SubscriptBox[\(u\), \(l\)] \(\*SubscriptBox[\(ϕ\), \(l\)]\)[x])\), \(2\)] \((\*UnderoverscriptBox[\(∑\), \(i = 1\), \(n\)]\*SubscriptBox[\(u\), \(i\)] \*SuperscriptBox[SubscriptBox[\(ω\), \(i\)], \(2\)] \(\*SubscriptBox[\(ϕ\), \(i\)]\)[x])\) \[DifferentialD]x\)\) - 
 Subscript[α, 2] Subscript[V, DC]^2 \!\(\*SubsuperscriptBox[\(∫\), \(0\), \(1\)]\(\(\*SubscriptBox[\(ϕ\), \(j\)]\)[x] \[DifferentialD]x\)\)) /.   j -> 1 ;

SysEq2 = (\!\(\*SubsuperscriptBox[\(∫\), \(0\), \(1\)]\(\(\*SubscriptBox[\(ϕ\), \(j\)]\)[x] \*SuperscriptBox[\((1 - \*UnderoverscriptBox[\(∑\), \(l = 1\), \(n\)]\*SubscriptBox[\(u\), \(l\)] \(\*SubscriptBox[\(ϕ\), \(l\)]\)[x])\), \(2\)] \((\*UnderoverscriptBox[\(∑\), \(i = 1\), \(n\)]\*SubscriptBox[\(u\), \(i\)] \*SuperscriptBox[SubscriptBox[\(ω\), \(i\)], \(2\)] \(\*SubscriptBox[\(ϕ\), \(i\)]\)[x])\) \[DifferentialD]x\)\) - 
 Subscript[α, 2] Subscript[V, DC]^2 \!\(\*SubsuperscriptBox[\(∫\), \(0\), \(1\)]\(\(\*SubscriptBox[\(ϕ\), \(j\)]\)[x] \[DifferentialD]x\)\)) /.   j -> 2 ;

SysEq3 = (\!\(\*SubsuperscriptBox[\(∫\), \(0\), \(1\)]\(\(\*SubscriptBox[\(ϕ\), \(j\)]\)[x] \*SuperscriptBox[\((1 - \*UnderoverscriptBox[\(∑\), \(l = 1\), \(n\)]\*SubscriptBox[\(u\), \(l\)] \(\*SubscriptBox[\(ϕ\), \(l\)]\)[x])\), \(2\)] \((\*UnderoverscriptBox[\(∑\), \(i = 1\), \(n\)]\*SubscriptBox[\(u\), \(i\)] \*SuperscriptBox[SubscriptBox[\(ω\), \(i\)], \(2\)] \(\*SubscriptBox[\(ϕ\), \(i\)]\)[x])\) \[DifferentialD]x\)\) - 
 Subscript[α, 2] Subscript[V, DC]^2 \!\(\*SubsuperscriptBox[\(∫\), \(0\), \(1\)]\(\(\*SubscriptBox[\(ϕ\), \(j\)]\)[x] \[DifferentialD]x\)\)) /.   j -> 3 ;

For[k = 2, k <= 100, k++,
Subscript[V, DC] = Subscript[V, dc][[k]]

Sol = Solve[{SysEq1 == 0, SysEq2 == 0, SysEq3 == 0}, {Subscript[u, 1],
 Subscript[u, 2], Subscript[u, 3]} ];

W = 0;
For[m = 1, m <= n, m++,
W = W + Subscript[u, m] *Subscript[ϕ, m][1] /. Sol[[1, m]];
];

W = Re[W ];
AppendTo[Solut, {Subscript[V, DC], W}];
];
Subscript[graphic, 1] =  ListPlot[Solut, PlotStyle -> PointSize[0.02], PlotRange -> All,  AxesLabel -> {"\!\(\*SubscriptBox[\(V\), \(DC\)]\)", 
"\!\(\*SubscriptBox[\(W\), \(Max\)]\)"}]

Answer 1

In this code we can replace Integrate with sufficient numeric integration using Gauss quadrature and Solve with NSolve, then we have:

Get["NumericalDifferentialEquationAnalysis`"];
np = 10; points = weights = Table[Null, {np}]; Do[
 points[[i]] = GaussianQuadratureWeights[np, 0, 1][[i, 1]], {i, 1, 
  np}];
Do[weights[[i]] = GaussianQuadratureWeights[np, 0, 1][[i, 2]], {i, 1, 
   np}];
GaussInt[f_, z_] := 
  Sum[(f /. z -> points[[i]])*weights[[i]], {i, 1, np}];
(*Geometrical data*)
d = 2*10^-6;  h = 2*10^-6;  b = 2*10^-6; l = 200*10^-6; EE = 
 1.66*10^11;  Eps = 8.854*10^-12;  ρ = 2332;

(*Numerical problem data*)
Subscript[α, 2] = (6*Eps*l^4)/(EE*h^3*
    d^3); Subscript[V, pull] = Sqrt[1.72/Subscript[α, 2]]; 
Subscript[V, dc] = (Range[100] - 1)*Subscript[V, pull]/100;
Solut = {};

(*Number of modes*)
n = 3; 

(*Normalized frequencies*)
Subscript[ω, 1] = 3.51602; Subscript[ω, 2] = 22.0345; 
Subscript[ω, 3] = 61.6972; Subscript[σ, 1] = 0.7341; 
Subscript[σ, 2] = 1.0185; Subscript[σ, 3] = 0.9992;

(*Modeshapes*)
Subscript[ϕ, j_][x_] := 
  Cosh[Sqrt[Subscript[ω, j]] x] - 
   Cos[Sqrt[Subscript[ω, j]] x] - 
   Subscript[σ, 
     j]*(Sinh[Sqrt[Subscript[ω, j]] x] -  
      Sin[Sqrt[Subscript[ω, j]] x]);

(*The derived eq for Cantilver beam *)
SysEq[j_] := (GaussInt[Subscript[ϕ, j][x] (1 - \!\(
\*UnderoverscriptBox[\(∑\), \(l = 1\), \(n\)]\(
\*SubscriptBox[\(u\), \(l\)] 
\(\*SubscriptBox[\(ϕ\), \(l\)]\)[x]\)\))^2 (\!\(
\*UnderoverscriptBox[\(∑\), \(i = 1\), \(n\)]\(
\*SubscriptBox[\(u\), \(i\)] 
\*SuperscriptBox[
SubscriptBox[\(ω\), \(i\)], \(2\)] 
\(\*SubscriptBox[\(ϕ\), \(i\)]\)[x]\)\)), x ] - 
     Subscript[α, 2] Subscript[V, DC]^2 GaussInt[
      Subscript[ϕ, j][x], x]) ;

For[k = 2, k <= 100, k++,
  Subscript[V, DC] = Subscript[V, dc][[k]];

  Sol = NSolve[Table[SysEq[j] == 0, {j, 3}], {Subscript[u, 1],
      Subscript[u, 2], Subscript[u, 3]} ];

  W = 0;
  For[m = 1, m <= n, m++,
   W = W + Subscript[u, m] *Subscript[ϕ, m][1] /. Sol[[1, m]];
   ];

  W = Re[W ];
  AppendTo[Solut, {Subscript[V, DC], W}];
  ];
Subscript[graphic, 1] =  
 ListPlot[Solut, PlotStyle -> PointSize[0.01], PlotRange -> All,  
  AxesLabel -> {"\!\(\*SubscriptBox[\(V\), \(DC\)]\)", 
    "\!\(\*SubscriptBox[\(W\), \(Max\)]\)"}]