Is this method allowed ?

$$left[begin{array}{rr|rr}
3 & 2 & 36 \
5 & 4 & 64
end{array}right]
sim
left[begin{array}{rr|rr}
1 & frac{2}{3} & 12 \
5 & 4 & 64
end{array}right]
sim left[begin{array}{rr|rr}
1 & frac{2}{3} & 12 \
0 & frac{2}{3} & 4
end{array}right] sim left[begin{array}{rr|rr}
1 & 0 & 8 \
0 & frac{2}{3} & 4
end{array}right] sim left[begin{array}{rr|rr}
1 & 0 & 8 \
0 & 1 & 6
end{array}right]
$$

which yields $x=8$ and $y=6$

The first step is $R_1 to R_1 times frac{1}{3}$

The second step is $R_2 to R_2 - 5R_1$

The third step is $R_1 to R_1 -R_2$

The fourth step is $R_2 to R_2times frac{3}{2}$

Here $R_i$ denotes the $i$ -th row.

How about using Cramer's Rule? Define $Delta_x=left[begin{matrix}36 & 2 \ 64 & 4end{matrix}right]$, $Delta_y=left[begin{matrix}3 & 36\ 5 & 64end{matrix}right]$
and $Delta_0=left[begin{matrix}3 & 2\ 5 &4end{matrix}right]$.

Now computation is trivial as you have: $x=dfrac{detDelta_x}{detDelta_0}$ and $y=dfrac{detDelta_y}{detDelta_0}$.

Fixed Point Iteration

This is not efficient but it's another valid way to solve the system. Treat the system as a matrix equation and rearrange to get $begin{bmatrix} x\ yend{bmatrix}$ on the left hand side.

Define
$fbegin{bmatrix} x\ yend{bmatrix}=begin{bmatrix} (36-2y)/3 \ (64-5x)/4end{bmatrix}$

Start with an intial guess of $begin{bmatrix} x\ yend{bmatrix}=begin{bmatrix} 0\ 0end{bmatrix}$

The result is $fbegin{bmatrix} 0\ 0end{bmatrix}=begin{bmatrix} 12\ 16end{bmatrix}$

Now plug that back into f

The result is $fbegin{bmatrix} 12\ 6end{bmatrix}=begin{bmatrix} 4/3\ 1end{bmatrix}$

Keep plugging the result back in. After 100 iterations you have:

$begin{bmatrix} 7.9991\ 5.9993end{bmatrix}$

Here is a graph of the progression of the iteration:

By false position:

Assume $x=10,y=3$, which fulfills the first equation, and let $x=10+x',y=3+y'$. Now, after simplification

$$3x'+2y'=0,\5x'+4y'=2.$$

We easily eliminate $y'$ (using $4y'=-6x'$) and get

$$-x'=2.$$

Though this method is not essentially different from, say elimination, it can be useful for by-hand computation as it yields smaller terms.

Another method to solve simultaneous equations in two dimensions, is by plotting graphs of the equations on a cartesian plane, and finding the point of intersection.

Any method you can come up with will in the end amount to Cramer's rule, which gives explicit formulas for the solution. Except special cases, the solution of a system is unique, so that you will always be computing the ratio of those determinants.

Anyway, it turns out that by organizing the computation in certain ways, you can reduce the number of arithmetic operations to be performed. For $2times2$ systems,
the different variants make little difference in this respect. Things become more interesting for $ntimes n$ systems.

Direct application of Cramer is by far the worse, as it takes a number of operations proportional to $(n+1)!$, which is huge. Even for $3times3$ systems, it should be avoided. The best method to date is Gaussian elimination (you eliminate one unknown at a time by forming linear combinations of the equations and turn the system to a triangular form). The total workload is proportional to $n^3$ operations.

The steps of standard Gaussian elimination:

$$begin{cases}ax+by=c,\dx+ey=f.end{cases}$$

Subtract the first times $dfrac da$ from the second,

$$begin{cases}ax+by=c,\0x+left(e-bdfrac daright)y=f-cdfrac da.end{cases}$$

Solve for $y$,

$$begin{cases}ax+by=c,\y=dfrac{f-cdfrac da}{e-bdfrac da}.end{cases}$$

Solve for $x$,

$$begin{cases}x=dfrac{c-bdfrac{f-cdfrac da}{e-bdfrac da}}a,\y=dfrac{f-cdfrac da}{e-bdfrac da}.end{cases}$$

So written, the formulas are a little scary, but when you use intermediate variables, the complexity vanishes:

$$d'=frac da,e'=e-bd',f'=f-cd'to y=frac{f'}{e'}, x=frac{c-by}a.$$

Anyway, for a $2times2$ system, this is worse than Cramer !

$$begin{cases}x=dfrac{ce-bf}{Delta},\y=dfrac{af-cd}{Delta}end{cases}$$ where $Delta=ae-bd$.

For large systems, say $100times100$ and up, very different methods are used. They work by computing approximate solutions and improving them iteratively until the inaccuracy becomes acceptable. Quite often such systems are sparse (many coefficients are zero), and this is exploited to reduce the number of operations. (The direct methods are inappropriate as they will break the sparseness property.)

Construct the Groebner basis of your system, with the variables ordered $x$, $y$:
$$ mathrm{GB}({3x+2y-36, 5x+4y-64}) = {y-6, x-8} $$
and read out the solution. (If we reverse the variable order, we get the same basis, but in reversed order.) Under the hood, this is performing Gaussian elimination for this problem. However, Groebner bases are not restricted to linear systems, so can be used to construct solution sets for systems of polynomials in several variables.

Perform lattice reduction on the lattice generated by $(3,2,-36)$ and $(5,4,-64)$. A sequence of reductions (similar to the Euclidean algorithm for GCDs): begin{align*}
(5,4,-64) - (3,2,-36) &= (2,2,-28) \
(3,2,-36) - (2,2,-28) &= (1,0,-8) tag{1} \
(2,2,-28) - 2(1,0,-8) &= (0,2,-12) tag{2} \
end{align*}
From (1), we have $x=8$. From (2), $2y = 12$, so $y = 6$. (Generally, there can be quite a bit more "creativity" required to get the needed zeroes in the lattice vector components. One implementation of the LLL algorithm, terminates with the shorter vectors ${(-1,2,4), (-2,2,4)}$, but we would continue to manipulate these to get the desired zeroes.)

$$begin{align}3x+2y&=36 tag1\ 5x+4y&=64tag2end{align}$$

From $(1)$, $x=frac{36-2y}{3}$, substitute in $(2)$ and you'll get $5(frac{36-2y}{3})+4y=64 implies y=6$ and then you can get that $x=24/3=8$

Another Method
From $(1)$, $x=frac{36-2y}{3}$

From $(2)$, $x=frac{64-4y}{5}$

But $x=x implies frac{36-2y}{3}=frac{64-4y}{5}$ do cross multiplication and you'll get $5(36-2y)=3(64-4y) implies y=6$ and substitute to get $x=8$

Other answers have given standard, elementary methods of solving simultaneous equations. Here are a few other ones that can be more long-winded and excessive, but work nonetheless.

Method $1$: (multiplicity of $y$)

Let $y=kx$ for some $kinBbb R$. Then $$3x+2y=36implies x(2k+3)=36implies x=frac{36}{2k+3}\5x+4y=64implies x(4k+5)=64implies x=frac{64}{4k+5}$$ so $$36(4k+5)=64(2k+3)implies (144-128)k=(192-180)implies k=frac34.$$ Now $$x=frac{64}{4k+5}=frac{64}{4cdotfrac34+5}=8implies y=kx=frac34cdot8=6.quadsquare$$

Method $2$: (use this if you really like quadratic equations :P)

How about we try squaring the equations? We get $$3x+2y=36implies 9x^2+12xy+4y^2=1296\5x+4y=64implies 25x^2+40xy+16y^2=4096$$ Multiplying the first equation by $10$ and the second by $3$ yields $$90x^2+120xy+40y^2=12960\75x^2+120xy+48y^2=12288$$ and subtracting gives us $$15x^2-8y^2=672$$ which is a hyperbola. Notice that subtracting the two linear equations gives you $2x+2y=28implies y=14-x$ so you have the nice quadratic $$15x^2-8(14-x)^2=672.$$ Enjoy!

As another iterative method I suggest the Jacobi Method. A sufficient criterion for its convergence is that the matrix must be diagonally dominant. Which this one in our system is not:

$begin{bmatrix} 3 &2 \ 5 &4end{bmatrix}begin{bmatrix} x \ yend{bmatrix}=begin{bmatrix}36 \ 64end{bmatrix}$

We can however fix this by replacing e.g. $y' := frac{1}{1.3} y$. Then the system is

$underbrace{begin{bmatrix} 3 & 2.6 \ 5 & 5.2end{bmatrix}}_{=:A}begin{bmatrix} x \ y'end{bmatrix}=begin{bmatrix}36 \ 64end{bmatrix}$

and $A$ is diagonally dominant. Then we can decompose $A = L + D + U$ into $L,U,D$ where $L,U$ are the strict upper and lower triangular parts and $D$ is the diagonal of $A$ and the iteration

$$vec x_{i+1} = - D^{-1}((L+R)vec x_i + b)$$

will converge to the solution $(x,y')$. Note that $D^{-1}$ is particularly easy to compute as you just have to invert the entries. So in theis case the iteration is

$$vec x_{i+1} = -begin{bmatrix} 1/3 & 0 \ 0 & 1/5.2 end{bmatrix}left(begin{bmatrix} 0 & 2.6 \ 5 & 0 end{bmatrix} vec x_i + bright)$$

So you can actually view this as a fixed point iteration of the function $f(vec x) = -D^{-1}((L+R)vec x + b)$ which is guaranteed to be a contraction in the case of diagonal dominance of $A$. It is actually quite slow and doesn't any practical application for directly solving systems of linear equations but it (or variations of it) is quite often used as a preconditioner.

It is clear that:

• 
  $x=10$, $y=3$ is an integer solution of $(1)$.


• 
  $x=12$, $y=1$ is an integer solution of $(2)$.

Then, from the theory of Linear Diophantine equations:

• Any integer solution of $(1)$ has the form $x_1=10+2t$, $y_1=3-3t$ with $t$ integer.


• Any integer solution of $(2)$ has the form $x_2=12+4t$, $y_2=1-5t$ with $t$ integer.

Then, the system has an integer solution $(x_0,y_0)$ if and only if there exists an integer $t$ such that

$$10+2t=x_0=12+4tqquadtext{and}qquad 3-3t=y_0=1-5t.$$

Solving for $t$ we see that there exists an integer $t$ satisfying both equations, which is $t=-1$. Thus the system has the integer solution
$$x_0=12+4(-1)=8,; y_0=1-5(-1)=6.$$

Note that we can pick any pair of integer solutions to start with. And the method will give the solution provided that the solution is integer, which is often not the case.

Consider the three vectors $textbf{A}=(3,2)$, $textbf{B}=(5,4)$ and $textbf{X}=(x,y)$. Your system could be written as $$textbf{A}cdottextbf{X}=a\textbf{B}cdottextbf{X}=b$$ where $a=36$, $b=64$ and $textbf{A}_{perp}=(-2,3)$ is orthogonal to $textbf{A}$. The first equation gives us $textbf{X}=dfrac{atextbf{A}}{textbf{A}^2}+lambdatextbf{A}_{perp}$. Now to find $lambda$ we use the second equation, we get $lambda=dfrac{b}{textbf{A}_{perp}cdottextbf{B}}-dfrac{atextbf{A}cdottextbf{B}}{textbf{A}^2timestextbf{A}_{perp}cdottextbf{B}}$. Et voilà :
$$textbf{X}=dfrac{atextbf{A}}{textbf{A}^2}+dfrac{textbf{A}_{perp}}{textbf{A}_{perp}cdottextbf{B}}left(b-dfrac{atextbf{A}cdottextbf{B}}{textbf{A}^2}right)$$