Solving nonlinear equation resulting from finite element method Using the finite element method (for a uniform mesh in the spatial domain) I have the system with initial conditions $u_j(0)=\cos(x_j)$ for $j=1,\dots, N$
$$\frac{d\vec{u}}{dt}=A\vec{u}+B\vec{c},$$
where $A$ and $B$ are coefficient matrices and $\vec{c}$ takes the form 
$$\vec{c}=\begin{bmatrix}\displaystyle\int_0^1\left(\vec{u}\cdot\beta(x)\right)^2\beta_j'(x)dx\end{bmatrix}_{1\leq j\leq N}$$
and $u_j=u_j(t)\approx u(x_j,t)$ is the semi-approximation of the solution at time $t$ on the mesh node $x_j$ for $j=1,\dots, N$ and $\beta_j$ is a linear basis function. My book says that this system can be solved with a "standard ODE solver", but provides no other insight besides that. My question is how does one go about solving such a system? Is there a particular MATLAB command that is used to solve this system? After searching online I haven't been able to find something conrete. Any discussion or insight would be greatly appreciated. 
 A: After doing the FEM discretization in space, you obtain a semi-discretized problem: Your temporal derivative $\frac{\mathrm d u }{\mathrm dt}$ remained, but you managed to approximate the spatial derivatives. The important thing now to realize is that your coefficients of the FEM functions are functions of time, i.e., $\vec{u} = \vec{u}(t)$.
To get a standard non-linear ODE, you want to get rid of the integral. You could do this in the standard procedure by transforming it to the reference element and then use a quadrature rule to transform the integral essentially into a sum.
You could also simplify the integral a bit a priori by noting that
$$\big(\vec{u}\cdot\beta(x)\big)^2 = \Bigg(\sum_k^N u_k \beta_k(x) \Bigg) \cdot\Bigg(\sum_i^N u_i \beta_i(x) \Bigg) = \sum_{i,k}^N u_ku_i\beta_k(x)\beta_i(x)$$
and thus
$$c_j = \int_0^1 \sum_{i,k}^N \beta_j'(x) u_ku_i\beta_k(x)\beta_i(x) \mathrm d x = \sum_{i,j}^N u_ju_i  \int_0^1 \beta_j'(x) \beta_k(x)\beta_i(x) \mathrm d x$$
This potentially allows you avoiding the quadrature since you know your basis functions $\vec{\beta}(x)$, i.e., you can precompute the integral.
In both cases, however, you end up with a problem like
$$\frac{\mathrm d \vec{u}}{\mathrm d t} = \vec{F}\big(\vec{u}(t)\big), $$
a nonlinear ODE for the coefficients $\vec{u}(t)$.
The standard command for this in Matlab would be ode45, although there are some other solvers for stiff systems, ...
