Solve equations with substitution method

Solve equations with substitution method is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Task

Let given equations:
3x + y = -1 and 2x - 3y = -19
Solve it with substitution method.

See related

• Reduced row echelon form
• Gaussian elimination

ALGOL 68

Translation of: Phix – second version (translation of Raku)

BEGIN # solve equations using the substitution method - translation of Phix  #

    PROC solve2 = ( []REAL e1, e2 )[]REAL:
         BEGIN
             REAL a1 = e1[ 1 ], b1 = e1[ 2 ], c1 = e1[ 3 ];
             REAL a2 = e2[ 1 ], b2 = e2[ 2 ], c2 = e2[ 3 ];
             REAL x  = (b2*c1 - b1*c2) / (b2*a1 - b1*a2);
             REAL y  = (a1*x - c1)/-b1;
             ( x, y )
         END # solve2 # ;

   OP   FMT = ( REAL v )STRING:     # prints v with at most 3 decimal places #
         BEGIN
            STRING result := fixed( ABS v, 0, 3 );
            IF result[ LWB result ] = "." THEN "0" +=: result FI;
            WHILE result[ UPB result ] = "0" DO result := result[ : UPB result - 1 ] OD;
            IF result[ UPB result ] = "." THEN result := result[ : UPB result - 1 ] FI;
            IF v < 0 THEN "-" ELSE "" FI + result
         END # FMT # ;

    []REAL xy = solve2( ( 3, 1, -1 ), ( 2, -3, -19 ) );
    print( ( "x = ", FMT xy[ 1 ], ", y = ", FMT xy[ 2 ], newline ) )

END

Output:

x = -2, y = 5

AWK

# syntax: GAWK -f SOLVE_EQUATIONS_WITH_SUBSTITUTION_METHOD.AWK
BEGIN {
    main("3,1,-1","2,-3,-19")
    exit(0)
}
function main(s1,s2,  arr,e1,e2,result_x,result_y,r1,r2,x1,x2,y1,y2) {
    split(s1,e1,",")
    split(s2,e2,",")
    x1 = e1[1]
    y1 = e1[2]
    r1 = e1[3]
    x2 = e2[1]
    y2 = e2[2]
    r2 = e2[3]
    arr[1] = x1
    arr[2] = -y1
    arr[3] = r1
    result_y = ((arr[1]*r2) - (x2*arr[3])) / ((x2*arr[2]) + (arr[1]*y2))
    result_x = (r1 - (y1*result_y)) / x1
    printf("x = %g\ny = %g\n",result_x,result_y)
}

Output:

x = -2
y = 5

BASIC

BASIC256

Translation of: FreeBASIC

arraybase 1
dim firstEquation(3)
firstEquation[1] = 3
firstEquation[2] = 1
firstEquation[3] = -1
dim secondEquation(3)
secondEquation[1] = 2
secondEquation[2] = -3
secondEquation[3] = -19

subroutine getCrossingPoint(firstEquation, secondEquation)
	x1 = firstEquation[1]
	y1 = firstEquation[2]
	r1 = firstEquation[3]
	x2 = secondEquation[1]
	y2 = secondEquation[2]
	r2 = secondEquation[3]
	dim temp(3)
	temp[1] =  x1
	temp[2] = -y1
	temp[3] =  r1
	resultY = ((temp[1]*r2) - (x2*temp[3])) / ((x2*temp[2]) + (temp[1]*y2))
	resultX = (r1 - (y1*resultY)) / x1
	print "x = "; resultX
	print "y = "; resultY
end subroutine

call getCrossingPoint(firstEquation, secondEquation)
end

Output:

Igual que la entrada de FreeBASIC.

FreeBASIC

Dim Shared As Integer firstEquation(1 To 3)  = { 3, 1, -1}
Dim Shared As Integer secondEquation(1 To 3) = { 2,-3,-19}

Sub getCrossingPoint(firstEquation() As Integer, secondEquation() As Integer)
    Dim As Integer x1 = firstEquation(1)
    Dim As Integer y1 = firstEquation(2)
    Dim As Integer r1 = firstEquation(3)
    Dim As Integer x2 = secondEquation(1)
    Dim As Integer y2 = secondEquation(2)
    Dim As Integer r2 = secondEquation(3)
    Dim As Integer temp(3)
    temp(1) =  x1
    temp(2) = -y1
    temp(3) =  r1
    Dim As Integer resultY = ((temp(1)* r2) - (x2 * temp(3))) / ((x2 * temp(2)) + (temp(1)*y2)) 
    Dim As Integer resultX = (r1 - (y1*resultY)) / x1 
    Print "x = "; resultX
	Print "y = "; resultY
End Sub

getCrossingPoint(firstEquation(), secondEquation())
Sleep

Output:

x = -2
y =  5

QBasic

Works with: QBasic

Works with: QuickBasic

Translation of: FreeBASIC

DIM firstEquation(3)
firstEquation(1) = 3
firstEquation(2) = 1
firstEquation(3) = -1
DIM secondEquation(3)
secondEquation(1) = 2
secondEquation(2) = -3
secondEquation(3) = -19

CALL getCrossingPoint(firstEquation(), secondEquation())
END

SUB getCrossingPoint (firstEquation(), secondEquation())
    x1 = firstEquation(1)
    y1 = firstEquation(2)
    r1 = firstEquation(3)
    x2 = secondEquation(1)
    y2 = secondEquation(2)
    r2 = secondEquation(3)    
    DIM temp(3)
    temp(1) = x1
    temp(2) = -y1
    temp(3) = r1
    resultY = ((temp(1) * r2) - (x2 * temp(3))) / ((x2 * temp(2)) + (temp(1) * y2))
    resultX = (r1 - (y1 * resultY)) / x1
    PRINT "x = "; resultX
    PRINT "y = "; resultY
END SUB

Output:

Igual que la entrada de FreeBASIC.

True BASIC

Works with: QBasic

Translation of: QBasic

DIM firstequation(3)
LET firstequation(1) = 3
LET firstequation(2) = 1
LET firstequation(3) = -1
DIM secondequation(3)
LET secondequation(1) = 2
LET secondequation(2) = -3
LET secondequation(3) = -19

SUB getcrossingpoint (firstequation(),secondequation())
    LET x1 = firstequation(1)
    LET y1 = firstequation(2)
    LET r1 = firstequation(3)
    LET x2 = secondequation(1)
    LET y2 = secondequation(2)
    LET r2 = secondequation(3)

    DIM temp(3)
    LET temp(1) = x1
    LET temp(2) = -y1
    LET temp(3) = r1

    LET resulty = ((temp(1)*r2)-(x2*temp(3)))/((x2*temp(2))+(temp(1)*y2))
    LET resultx = (r1-(y1*resulty))/x1
    PRINT "x = "; resultx
    PRINT "y = "; resulty
END SUB

CALL getcrossingpoint (firstequation(), secondequation())
END

Output:

Igual que la entrada de FreeBASIC.

Yabasic

Translation of: FreeBASIC

dim firstEquation(3)
firstEquation(1) = 3
firstEquation(2) = 1
firstEquation(3) = -1
dim secondEquation(3)
secondEquation(1) = 2
secondEquation(2) = -3
secondEquation(3) = -19

sub getCrossingPoint(firstEquation(), secondEquation())
    x1 = firstEquation(1)
    y1 = firstEquation(2)
    r1 = firstEquation(3)
    x2 = secondEquation(1)
    y2 = secondEquation(2)
    r2 = secondEquation(3)
    dim temp(3)
    temp(1) =  x1
    temp(2) = -y1
    temp(3) =  r1
    resultY = ((temp(1)*r2) - (x2*temp(3))) / ((x2*temp(2)) + (temp(1)*y2)) 
    resultX = (r1 - (y1*resultY)) / x1 
    print "x = ", resultX
    print "y = ", resultY
end sub

getCrossingPoint(firstEquation(), secondEquation())
end

Output:

Igual que la entrada de FreeBASIC.

Delphi

Works with: Delphi version 6.0

Library: SysUtils,StdCtrls

Rather than just provide a one-off solution for this one problem, I've included code that will solve any system of equations with n-equations and n-unknowns. The code revolves around a matrix object, that contains an N by N matrix. The object can perform Gausian Elimination, which reduces the matrix to "Row Eschelon" format. From Row Eschelon format, you can back substitute and get the values for the unknowns. It also includes code to do "Gauss-Jordan" elimination, which eliminates the back substitution step. The object presented here is a subset of a more elaborate object that can perform all kinds of matrix operations.

{This code normally resides in a library, it is provided here for clarity}

type T2DMatrix = class(TObject)
 private
  FDoubleArray: array of double;
  FRows,FColumns: integer;
  procedure Put(Row,Col: integer; Item: double);
  function Get(Row,Col: integer): double;
 protected
 public
  constructor Create(Row,Col: integer);
  property Rows: integer read FRows;
  property Columns: integer read FColumns;
  procedure SetDimensions(Row,Col: integer);
  property Value[Row,Col: integer]: Double read Get write Put; default;
  function BackSubstitute: TDoubleDynArray;
  procedure GausianElimination;
  function GaussBackSubsituation: TDoubleDynArray;
  procedure GaussJordanElimination;
  procedure ExchangeRows(R1, R2: integer);
  function MatrixStr(Digits: integer): string;
 end;


{ T2DMatrix }


procedure T2DMatrix.SetDimensions(Row, Col: integer);
begin
FRows:=Row; FColumns:=Col;
SetLength(FDoubleArray,Row * Col);
end;

constructor T2DMatrix.Create(Row, Col: integer);
begin
SetDimensions(Row, Col);
end;


procedure T2DMatrix.Put(Row,Col: integer; Item: double);
{Insert Double at index}
var Inx: integer;
begin
Inx:=(Col * Rows) + Row;
FDoubleArray[Inx]:=Item;
end;


function T2DMatrix.Get(Row,Col: integer): double;
{Get Double at the index}
var Inx: integer;
begin
Inx:=(Col * Rows) + Row;
Result:=FDoubleArray[Inx];
end;

{Matrix operations}

procedure T2DMatrix.ExchangeRows(R1,R2: integer);
{Exchange the specified Rows}
var Col: integer;
var T: double;
begin
for Col:=0 to Self.Columns-1 do
	begin
	T:=Self[R1,Col];
	Self[R1,Col]:=Self[R2,Col];
	Self[R2,Col]:=T;
	end;
end;



procedure T2DMatrix.GausianElimination;
var I,K,J : Integer;
var S : double;
{Small value to prevent divide by zero}
const Epsilon = 5.0e-162;
begin
{Do gaussian elimination and convert Row Echelon}
K:=0;
while True do
	begin
	{If the pivot is zero, find another row with non-zero in the same column}
	if Self[K,K]=0 then
		begin
		for I:=K+1 to Self.Rows-1 do
		 if Self[I,K]<>0 then
		 	begin
		 	Self.ExchangeRows(K,I);
		 	break;
		 	end;
		 end;

	{Get pivot again}
	S:=Self[K,K];
	{if it is still zero, prevent divide by zero}
	if S=0.0 then S:=Epsilon;
	{Use "scaling primative row operation" to set pivot to one i.e divide each item by pivot}
	for J:=K to Self.Columns-1 do Self[K,J]:=Self[K,J]/S;
	{Exit if we are on the bottom row}
	if K>=Self.Rows-1 then break;
	{Now that the previous row has 1 in the leading coefficient (Pivot)}
	{We convert all remaining columns below this item to zero }
	{Using "pivot primative row operation" = Current Row - Pivot * Start Row}
	for I:=K+1 to Self.Rows-1 do
		{Do this to current column for all remaining rows}
		begin
		{Get the leading coefficient}
		S:=Self[I,K];
		{Use it to zero column and apply to all items in row}
		for J:=K to Self.Columns-1 do Self[I,J]:=Self[I,J]-S*Self[K,J];
		end;
	{Point to next pivot}
	K:=K+1;
	end;
{Matrix is now in Row Echelon form }
end;




function T2DMatrix.BackSubstitute: TDoubleDynArray;
var Row,J: integer;
var Sum: double;
begin
SetLength(Result,Self.Rows);
for Row:=Self.Rows-1 downto 0 do
	begin
	Result[Row]:=Self[Row,Self.Columns-1];
	Sum:=0;
	for J:=Self.Rows-1 downto Row+1 do
        Sum:=Sum + Result[J] * Self[Row,J];
	Result[Row]:=Result[Row] - Sum;
	end;
end;



function T2DMatrix.GaussBackSubsituation: TDoubleDynArray;
begin
Self.GausianElimination;
Result:=Self.BackSubstitute;
end;



procedure T2DMatrix.GaussJordanElimination;
{ Do Gauss Jordan Elimination }
var I,K,J : Integer;
var S : double;
var X,Y: integer;
var Str: string;
{Small value to prevent divide by zero}
const Epsilon = 5.0e-162;
begin
{Do Gaussian Elimination to put matrix in Row Echelon format}
Self.GausianElimination;
{Now do the Jordan part to put it in Reduced Row Echelon form}
{Point K at the both the source row and target column (left part of matrix is square)}
K:=Self.Rows-1;
while K>0 do
	begin
	{Point I at the target row}
	I:=K-1;
	repeat
		begin
		{Get the coefficient of target column and row, the value we want to zero }
		S:=Self[I,K];
		{Multiply source row by the target coefficient and subtract from each item in the target row}
		{Since the target column is one in the source row, this will zero out the target coefficient}
		for J:=K to Self.Columns-1 do Self[I,J]:=Self[I,J]-S*Self[K,J];
		{Point to previous row}
		I:=I-1;
		end
	until I<0;
	{Point to previous column}
	K:=K-1;
	end;
end;



function T2DMatrix.MatrixStr(Digits: integer): string;
var Row,Col: integer;
begin
Result:=IntToStr(Self.Rows)+'X'+IntToStr(Self.Columns)+CRLF_Char;
for Row:=0 to Self.Rows-1 do
	begin
	Result:=Result+'[';
	for Col:=0 to Self.Columns-1 do
		begin
		if Col<>0 then Result:=Result+'  ';
		Result:=Result+FloatToStrF(Self[Row,Col],ffFixed,18,Digits);
		end;
	Result:=Result+']'+CRLF_Char;
	end;
end;

{===============================================================================}

procedure SolveEquations(Memo: TMemo);
var Mat: T2DMatrix;
var DA: TDoubleDynArray;
var I: integer;
var S: string;
begin
Mat:=T2DMatrix.Create(2,3);
try
{3x + y = -1}
{2x - 3y = -19}
Mat[0,0]:=3; Mat[0,1]:=1; Mat[0,2]:=-1;
Mat[1,0]:=2; Mat[1,1]:=-3; Mat[1,2]:=-19;
Memo.Lines.Add('Solve with Gaussian Elimination and Substitution');
Memo.Lines.Add('------------------------------------------------');
Memo.Lines.Add(Mat.MatrixStr(1));
Mat.GausianElimination;
Memo.Lines.Add('Row Echelon after Gaussian Elimination');
Memo.Lines.Add(Mat.MatrixStr(1));
DA:=Mat.BackSubstitute;
Memo.Lines.Add('Matrix after Back Substitution');
Memo.Lines.Add(Mat.MatrixStr(1));
Memo.Lines.Add(Format('Solution: X=%2.2f, Y=%2.2f',[DA[0],DA[1]]));

Mat[0,0]:=3; Mat[0,1]:=1; Mat[0,2]:=-1;
Mat[1,0]:=2; Mat[1,1]:=-3; Mat[1,2]:=-19;

Memo.Lines.Add('');
Memo.Lines.Add('Solve with Gaussian Jordan elimination');
Memo.Lines.Add('--------------------------------------');
Memo.Lines.Add(Mat.MatrixStr(1));
Mat.GaussJordanElimination;
Memo.Lines.Add('Matrix after Gauss-Jordan');
Memo.Lines.Add(Mat.MatrixStr(1));
Memo.Lines.Add(Format('Solution: X=%2.2f, Y=%2.2f',[Mat[0,2],Mat[1,2]]));
finally Mat.Free; end;
end;

Output:

Solve with Gaussian Elimination and Substitution
------------------------------------------------
2X3
[3.0  1.0  -1.0]
[2.0  -3.0  -19.0]

Row Echelon after Gaussian Elimination
2X3
[1.0  0.3  -0.3]
[0.0  1.0  5.0]

Matrix after Back Substitution
2X3
[1.0  0.3  -0.3]
[0.0  1.0  5.0]

Solution: X=-2.00, Y=5.00

Solve with Gaussian Jordan elimination
--------------------------------------
2X3
[3.0  1.0  -1.0]
[2.0  -3.0  -19.0]

Matrix after Gauss-Jordan
2X3
[1.0  0.0  -2.0]
[0.0  1.0  5.0]

Solution: X=-2.00, Y=5.00
Elapsed Time: 21.255 ms.

EasyLang

func[] solve a[] b[] .
   r[] &= (b[2] * a[3] - a[2] * b[3]) / (b[2] * a[1] - a[2] * b[1])
   r[] &= (a[1] * r[1] - a[3]) / -a[2]
   return r[]
.
print solve [ 3 1 -1 ] [ 2 -3 -19 ]

Output:

[ -2 5 ]

jq

Works with jq, the C implementation of jq

Works with gojq, the Go implementation of jq

The solution presented here handles all the edge cases.

With trivial modifications, the following will also work with jaq, the Rust implementation of jq.

# The equation ax + by = c is represented by the array [a, b, c] 
def solve( $e1; $e2 ):
    $e1 as [$a1, $b1, $c1]
  | $e2 as [$a2, $b2, $c2]
  | ($b2 * $a1  -  $b1 * $a2 ) as $d
  | if $d == 0 then "there is no unique solution as the discriminant is 0" | error
    else
      {  x : (($b2 * $c1  -  $b1 * $c2 ) / $d) }
      | if   $b1 != 0
        then .y = ( $a1 * .x   -  $c1 ) / -$b1
        else .y = ( $a2 * .x   -  $c2 ) / -$b2
        end
    end;

solve( [3,1,-1]; [2,-3,-19] )

Output:

{
  "x": -2,
  "y": 5
}

Julia

function parselinear(s)
    ab, c = strip.(split(s, "="))
    a, by = strip.(split(ab, "x"))
    b = replace(by, r"[\sy]" => "")
    b[end] in "-+" && (b *= "1")
    b = replace(b, "--" => "")
    return map(x -> parse(Float64, x == "" ? "1" : x), [a, b, c])
end

function solvetwolinear(s1, s2)
    a1, b1, c1 = parselinear(s1)
    a2, b2, c2 = parselinear(s2)
    x = (b2 * c1 - b1 * c2) / (b2 * a1 - b1 * a2)
    y = (a1 * x - c1 ) / -b1
    return x, y
end

@show solvetwolinear("3x + y = -1", "2x - 3y = -19")  # solvetwolinear("3x + y = -1", "2x - 3y = -19") = (-2.0, 5.0)

Nim

Translation of: Python

type Equation = tuple[cx, cy, cr: float]   # cx.x + cy.y = cr.

func getCrossingPoint(firstEquation, secondEquation: Equation): tuple[x, y: float] =
  let (x1, y1, r1) = firstEquation
  let (x2, y2, r2) = secondEquation
  let temp = (x1, -y1, r1)
  result.y = ((temp[0] * r2) - (x2 * temp[2])) / ((x2 * temp[1]) + (temp[0] * y2))
  result.x = (r1 - (y1 * result.y)) / x1

when isMainModule:
  let firstEquation: Equation = (3, 1, -1)
  let secondEquation: Equation = (2, -3, -19)
  let (x, y) = getCrossingPoint(firstEquation, secondEquation)
  echo "x = ", x
  echo "y = ", y

Output:

x = -2.0
y = 5.0

Perl

use strict;
use warnings;
use feature 'say';

sub parse {
    my($e) = @_;
    $e =~ s/ ([xy])/ 1$1/;
    $e =~ s/[ =\+]//g;
    split /[xy=]/, $e;
}

sub solve {
    my($a1, $b1, $c1, $a2, $b2, $c2) = @_;
    my $X = ( $b2 * $c1  -  $b1 * $c2 )
          / ( $b2 * $a1  -  $b1 * $a2 );
    my $Y = ( $a1 * $X  -  $c1 ) / -$b1;
    return $X, $Y;
}

say my $result = join ' ', solve( parse('3x + y = -1'), parse('2x - 3y = -19') );

Output:

-2 5

Phix

Slightly modified copy of solveN() from Solving_coin_problems#Phix, admittedly a tad overkill for this task, as it takes any number of rules and any number of variables.

with javascript_semantics
procedure solve(sequence rules, unknowns)
--
-- Based on https://mathcs.clarku.edu/~djoyce/ma105/simultaneous.html
--  aka the ancient Chinese Jiuzhang suanshu ~100 B.C. (!!)
--
-- Example:
--  rules = {{18,1,1},{38,1,5}}, ie 18==p+n, 38==p+5*n
--  unknowns = {"pennies","nickels"}
--
--  In the elimination phase, both p have multipliers of 1, so we can
--  ignore those two sq_mul and just do (38=p+5n)-(18=p+n)==>(20=4n).
--  Obviously therefore n is 5 and substituting backwards p is 13.
--
    string res
    sequence sentences = rules, ri, rj
    integer l = length(rules), rii, rji
    rules = deep_copy(rules)
    for i=1 to l do
        -- successively eliminate (grow lower left triangle of 0s)
        ri = rules[i]
        if length(ri)!=l+1 then ?9/0 end if
        rii = ri[i+1]
        if rii=0 then ?9/0 end if
        for j=i+1 to l do
            rj = rules[j]
            rji = rj[i+1]
            if rji!=0 then
                rj = sq_sub(sq_mul(rj,rii),sq_mul(ri,rji))
                if rj[i+1]!=0 then ?9/0 end if -- (job done)
                rules[j] = rj
            end if
        end for 
    end for 
    for i=l to 1 by -1 do
        -- then substitute each backwards
        ri = rules[i]
        rii = ri[1]/ri[i+1] -- (all else should be 0)
        rules[i] = sprintf("%s = %d",{unknowns[i],rii})
        for j=i-1 to 1 by -1 do
            rj = rules[j]
            rji = rj[i+1]
            if rji!=0 then
                rules[j] = 0
                rj[1] -= rji*rii
                rj[i+1] = 0
                rules[j] = rj
            end if
        end for
    end for 
    res = join(rules,", ")
    printf(1,"%v ==> %s\n",{sentences,res})
end procedure

--for 3x + y = -1 and 2x - 3y = -19:
solve({{-1,3,1},{-19,2,-3}},{"x","y"})

Output:

{{-1,3,1},{-19,2,-3}} ==> x = -2, y = 5

Alternatively, since I'm staring right at it, here's a

Translation of: Raku

with javascript_semantics
function solve2(sequence e1,e2)
    atom {a1,b1,c1} = e1,
         {a2,b2,c2} = e2,
         x = (b2*c1 - b1*c2)
           / (b2*a1 - b1*a2),
         y = (a1*x - c1)/-b1;
    return {x, y}
end function
printf(1,"x = %d, y = %d\n",solve2({3,1,-1},{2,-3,-19}))

Output:

x = -2, y = 5

Python

#!/usr/bin/python

firstEquation  = [ 3, 1, -1]
secondEquation = [ 2,-3,-19]

def getCrossingPoint(firstEquation, secondEquation):
    x1 = firstEquation[0]
    y1 = firstEquation[1]
    r1 = firstEquation[2]
    x2 = secondEquation[0]
    y2 = secondEquation[1]
    r2 = secondEquation[2]
    temp = []
    temp.append( x1)
    temp.append(-y1)
    temp.append( r1)
    resultY = ((temp[0]*r2) - (x2*temp[2])) / ((x2*temp[1]) + (temp[0]*y2)) 
    resultX = (r1 - (y1*resultY)) / x1 
    print("x = ", resultX)
    print("y = ", resultY)


if __name__ == "__main__":
    getCrossingPoint(firstEquation, secondEquation)

Output:

x =  -2.0
y =  5.0

Raku

sub solve-system-of-two-linear-equations ( [ \a1, \b1, \c1 ], [ \a2, \b2, \c2 ] ) {
    my \X = ( b2 * c1   -   b1 * c2 )
          / ( b2 * a1   -   b1 * a2 );

    my \Y = ( a1 * X    -   c1 ) / -b1;

    return X, Y;
}
say solve-system-of-two-linear-equations( (3,1,-1), (2,-3,-19) );

Output:

(-2 5)

Ring

firstEquation = [3.0,1.0,-1.0] secondEquation = [2.0,-3.0,-19.0]
getCrossingPoint(firstEquation,secondEquation)

func getCrossingPoint(firstEquation,secondEquation)
     x1 = firstEquation[1] y1 = firstEquation[2] r1 = firstEquation[3] x2 = secondEquation[1] y2 = secondEquation[2] r2 = secondEquation[3]
     temp = []
     add(temp,x1) add(temp,-y1) add(temp,r1)
     resultY = ((temp[1]* r2) - (x2 * temp[3])) / ((x2 * temp[2]) + (temp[1]*y2)) resultX = (r1 - (y1*resultY)) / x1 
     see "x = " + resultX + nl + "y = " + resultY + nl

Output:

x = -2
y = 5

Wren

var solve = Fn.new { |e1, e2|
    e2 = e2.toList
    for (i in 1..2) e2[i] = e2[i] * e1[0] / e2[0]
    var y = (e2[2] - e1[2]) / (e2[1] - e1[1])
    var x = (e1[2] - e1[1] * y) / e1[0]
    return [x, y]
}

var e1 = [3, 1, -1]
var e2 = [2, -3, -19]
var sol = solve.call(e1, e2)
System.print("x = %(sol[0]), y = %(sol[1])")

Output:

x = -2, y = 5

XPL0

This shows the vector routines from xpllib.xpl.

func real VSub(A, B, C);        \Subtract two 3D vectors
real    A, B, C;                \A:= B - C
[A(0):= B(0) - C(0);            \VSub(A, A, C) => A:= A - C
 A(1):= B(1) - C(1);
 A(2):= B(2) - C(2);
return A;
];      \VSub

func real VMul(A, B, S);        \Multiply 3D vector by a scalar
real    A, B, S;                \A:= B * S
[A(0):= B(0) * S;               \VMul(A, A, S) => A:= A * S
 A(1):= B(1) * S;
 A(2):= B(2) * S;
return A;
];      \VMul

real E1, E2, X1, X2, X, Y;
[E1:= [3.,  1.,  -1.];
 E2:= [2., -3., -19.];
X1:= E1(0);
X2:= E2(0);
VMul(E1, E1, X2);
VMul(E2, E2, X1);
VSub(E1, E1, E2);
Y:= E1(2)/E1(1);
E2(1):= E2(1)*Y;
E2(2):= E2(2)-E2(1);
X:= E2(2)/E2(0);
Text(0, "x = ");  RlOut(0, X);  CrLf(0);
Text(0, "y = ");  RlOut(0, Y);  CrLf(0);
]

Output:

x =    -2.00000
y =     5.00000