Problem description: I want to numerically solve system of two time-independent partial differential equations (pde) of the following simple form
$$\frac{\partial u(x,y)}{\partial x} = f_1(x,y,u),$$
$$\frac{\partial u(x,y)}{\partial y} = f_2(x,y,u),$$
where $u(x,y)$ is the unknown we want to solve for, $x$ and $y$ are the spatial variables and $f_1$,$f_2$ are some smooth (or at least continuously differentiable) and possibly nonlinear functions satisfying $\frac{\partial f_1}{\partial y}=\frac{\partial f_2}{\partial x}$. More compactly
$$\nabla u(x,y) = f(x,y,u),$$
with $\nabla:=[\frac{\partial~ \cdot}{\partial x}, \frac{\partial~ \cdot}{\partial y}]^T$ and $f:=[f_1,f_2]^T$.
I want to solve it on a rectangular domain (defined let’s say by $x$,$y$ for which $W>x>0$ and $H>y>0,~ W,H \in \mathbb R$) for an initial given boundary condition of $u(0,0)=c,~c \in \mathbb R$.
What I have tried so far: $u(x,y)$ can be imagined like a surface. We stand at the initial point and the partial derivatives tell us how much the surface goes upwards or downwards in the corresponding directions. In 1D case, when we have $\frac{du}{dx}=f(x,u)$, an ODE solver (like ode45 in Matlab) can be used (we can interpret the time $t$ as the spatial variable $x$ without no harm here). I thought, whether there is some way how to extend this approach to the 2 dimensional case. The solution that works to some extent is to solve for $u$ along one of the edges of the domain (e.g. find $u$ on $W>x>0,~y=0$ using ode45 as mentioned above) and then to use points of this solution as starting values for other set of ode's (now going in the perpendicular direction). In other words, if the goal is to find the values of $u$ on a regular grid, then we first find its values along the first row of the grid and then use the result as a starting point for ode's going along each of the columns and thus getting values for the whole grid. Although such a solution gives an idea how the surface $u(x,y)$ looks like, it is, however, not very precise as the individual solutions for rows/columns of the grid get apart from each other as we are moving away from the edge used as a set of initial values.
Note: I have also tried chebfun toolbox for Matlab and pde toolbox, but it seems that they only support time-dependent equations.