diff --git a/Mex/a21.cpp b/Mex/a21.cpp
index 905870c..4f3b314 100644
--- a/Mex/a21.cpp
+++ b/Mex/a21.cpp
@@ -46,11 +46,12 @@ void mexFunction(int nlhs, mxArray *plhs[], 			// Output variables
 				x[i] = x[i] - GA[i][j] * x[j];
 		x[i] = x[i] / GA[i][i];							//now finally divide the rhs by the coefficient of the variable to be calculated
 	}
-	
+	nlhs = variables;
 	for(int j=0; j < variables; j++){
 		plhs[j]=mxCreateDoubleScalar(x[j]);
 		
 	}
+	
 	for (int i = 0; i <= variables; i++ ){
 		free(GA[i]);
 	}
diff --git a/Mex/a40.cpp b/Mex/a40.cpp
new file mode 100644
index 0000000..4a5585d
--- /dev/null
+++ b/Mex/a40.cpp
@@ -0,0 +1,91 @@
+#include "mex.h" 
+#include "matrix.h"
+#include "stdlib.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], 			// Output variables 
+				 int nrhs, const mxArray *prhs[]) 		// Input variables 
+{
+	int variables = *mxGetPr(prhs[0]);
+	
+	double* x = (double *)mxCalloc(variables+1, sizeof(double));					//create solution array
+	
+	double** GA = (double **)mxCalloc(variables+1, sizeof(double*));				//create Gauss Array (GA)
+	for (int i = 0; i <= variables+1; i++){
+		GA[i] = (double *)mxCalloc(variables+1, sizeof(double));
+	}
+	
+	for(int i = 0; i < variables*(variables+1); i++){							//copy input array into GA								
+		GA[i%variables][i/variables] = mxGetPr(prhs[1])[i];	
+	}
+	
+	
+	int i, j, k;
+	for (i = 0; i<variables; i++)						//Pivotisation
+		for (k = i + 1; k<variables; k++)
+			if (abs(GA[i][i])<abs(GA[k][i]))
+				for (j = 0; j <= variables; j++)
+				{
+					long double temp = GA[i][j];
+					GA[i][j] = GA[k][j];
+					GA[k][j] = temp;
+				}
+
+	for (i = 0; i<variables - 1; i++)					//gauss elimination
+		for (k = i + 1; k<variables; k++)
+		{
+			long double t = GA[k][i] / GA[i][i];
+			for (j = 0; j <= variables; j++)
+				GA[k][j] = GA[k][j] - t*GA[i][j];		//make the elements below the pivot elements equal to zero or elimnate the variables
+		}
+
+	for (i = variables - 1; i >= 0; i--)                //back-substitution
+	{													//x is an array whose values correspond to the values of x,y,z..
+		x[i] = GA[i][variables];						//make the variable to be calculated equal to the rhs of the last equation
+		for (j = i + 1; j<variables; j++)
+			if (j != i)									//then subtract all the lhs values except the coefficient of the variable whose value is being calculated
+				x[i] = x[i] - GA[i][j] * x[j];
+		x[i] = x[i] / GA[i][i];							//now finally divide the rhs by the coefficient of the variable to be calculated
+	}
+	
+	nlhs = variables;
+	for(int l=0;l<variables;l++){
+		char temp;
+		
+		if(x[variables-1-l] != 0)
+		{
+			if(l!=0){
+			mexPrintf("+");
+			}
+			if(x[variables-1-l] < 0 ){
+				mexPrintf("-");
+				x[variables-l] *= -1;
+			}
+			temp=(char)(48+x[variables-1-l]);
+			mexPrintf("%c",temp);
+			
+			if(variables-1-l != 0){
+				mexPrintf("*");
+				mexPrintf("x");
+				mexPrintf("^");
+				temp=(char)(48+variables-1-l);
+				mexPrintf("%c",temp);
+			}
+			
+		}
+		
+	}
+	mexPrintf("\n");
+	for(int j=0; j < variables; j++){
+		plhs[j]=mxCreateDoubleScalar(x[j]);
+		
+	}
+	
+	for (int i = 0; i <= variables; i++ ){
+		free(GA[i]);
+	}
+	
+	free(GA);
+	free(x);
+  
+return;
+}
\ No newline at end of file
diff --git a/Mex/a40.m b/Mex/a40.m
new file mode 100644
index 0000000..432eb93
--- /dev/null
+++ b/Mex/a40.m
@@ -0,0 +1,18 @@
+% function accepts 1 double and 1 double matrix (#of_variables, [row1;row2])
+% each row1 contains #of_variables+1 elements
+% returns solution vector
+% outputs the function
+% bsp.:
+%			6x+12y=30
+%			3x+ 3y= 9
+%		a21(2,[6,12,30;3,3,9]);
+%			2*x^1+1				<-function
+%
+%				OR
+%
+%			6x+12y=30
+%			3x+ 3y= 9
+%		[a1,a2]=a21(2,[6,12,30;3,3,9])
+%			2*x^1+1				<-function
+%			a1 = 1				<-factors of x^...
+%			a2 = 2
\ No newline at end of file
diff --git a/Mex/a700.cpp b/Mex/a700.cpp
new file mode 100644
index 0000000..0cfcf35
--- /dev/null
+++ b/Mex/a700.cpp
@@ -0,0 +1,53 @@
+#include "mex.h" 
+#include "matrix.h"
+#include "stdlib.h"
+
+
+
+void mexFunction(int nlhs, mxArray *plhs[], 			// Output  
+				 int nrhs, const mxArray *prhs[]) 		// Input  
+{
+	char *func_Name;
+	func_Name = mxArrayToString(prhs[0]);
+	double t0 = *mxGetPr(prhs[1]);
+	double tfinal= *mxGetPr(prhs[2]);
+	double x0 = *mxGetPr(prhs[3]);
+	double* tout = mxGetPr(prhs[4]);
+	double* xout = mxGetPr(prhs[5]);
+	
+	double 	theta_umg = 20,								//temp of room
+			theta = x0,									//start temp
+			h = 1,										//incremtens size
+			t_start = t0, 								//start time
+			t_end = tfinal, 							//end time
+			t=t_start,									
+			theta_dot;
+	double* t_out = (double *)mxCalloc((t_end-t_start)/h, sizeof(double));			//array cont. time stamps
+	double* theta_out = (double *)mxCalloc((t_end-t_start)/h, sizeof(double));		//array cont. temp
+
+	t_out[0]=t_start;
+	theta_out[0]=theta;
+	
+	for(int i = 1; t<t_end; i++)
+	{
+		theta_dot = -1/2*(theta-theta_umg);
+		theta += theta_dot*h;
+		t += h;
+		t_out[i] = t;
+		theta_out = [i]= theta;
+	}
+	
+
+	mxArray *func_Input[2];
+	mxArray *func_Outputs[1]; 
+
+	func_Input[0] = ;
+	func_Input[1] = ;
+	mexCallMATLAB(1, func_Outputs , 1 , func_Input , "plot");
+
+// Output
+
+
+		 
+ return;
+}
\ No newline at end of file