diff --git a/lab3/omp_heat2D.c b/lab3/omp_heat2D.c
new file mode 100644
index 0000000..b4d9ed9
--- /dev/null
+++ b/lab3/omp_heat2D.c
@@ -0,0 +1,232 @@
+/****************************************************************************
+ * DESCRIPTION:  
+ *   Serial HEAT2D Example - C Version
+ *   This example is based on a simplified 
+ *   two-dimensional heat equation domain decomposition.  The initial 
+ *   temperature is computed to be high in the middle of the domain and 
+ *   zero at the boundaries.  The boundaries are held at zero throughout 
+ *   the simulation.  During the time-stepping, an array containing two 
+ *   domains is used; these domains alternate between old data and new data.
+ *
+ *  The physical region, and the boundary conditions, are suggested
+    by this diagram;
+
+                   u = 0
+             +------------------+
+             |                  |
+    u = 100  |                  | u = 100
+             |                  |
+             +------------------+
+                   u = 100
+
+Interrior point :
+  u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
+
+ ****************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <omp.h>
+
+#define N 500
+#define M 500  
+
+#define ITER_PRINT 100
+#define PRINT_DATA 1
+
+#define EPSILON 1e-1
+
+
+void update(int nx,int ny, float *u, float *unew, float * diff);
+void inidat(int nx, int ny, float *u, float *unew); 
+void prtdat(int nx, int ny, float *u,const char *fnam);
+
+
+
+
+int main(int argc, char *argv[])
+{
+
+    float diff=1.0;
+
+    float *u     = (float *)malloc(N * M * sizeof(float));
+    float *unew  = (float *)malloc(N * M * sizeof(float));
+ 
+    if(u==0 || unew ==0)
+    {
+        perror("Can't allocated data\n");
+        return -1;
+    }
+
+     printf ( "\n" );
+     printf ( "HEATED_PLATE\n" );
+     printf ( "  Parallel OpenMP version, using %d Threads\n",omp_get_max_threads() );
+     printf ( "  A program to solve for the steady state temperature distribution\n" );
+     printf ( "  over a rectangular plate.\n" );
+     printf ( "  Spatial grid of %d by %d points.\n\n", M, N );
+
+
+    /* Initialize grid and create input file */
+    printf("Initializing grid\n");
+    
+    inidat(N, M,u,unew);
+
+    prtdat(N, M,u, "initial.dat");
+    
+   
+    printf("Start computing\n");
+
+    int iter=0;
+
+    /* 
+     *   iterate until the  new solution unew differs from the old solution u
+     *     by no more than EPSILON.
+     *     */
+     
+    while(diff> EPSILON) {
+
+        update(N, M, u, unew,&diff);
+    
+        if(iter%ITER_PRINT==0)
+        printf("Iteration %d, diff = %f\n ", iter,diff);
+
+        iter++;
+    }
+
+    prtdat(N, M, u, "final.dat");
+     
+    free(u);
+    free(unew);
+}
+
+
+
+/****************************************************************************
+ *  subroutine update
+ ****************************************************************************/
+void update(int nx,int ny, float *u, float *unew, float * diff)
+{
+    int ix, iy;
+    *diff=0.0;
+
+#pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy)
+    for (ix = 1; ix < nx-1; ix++) {
+        for (iy = 1; iy < ny-1; iy++) {
+            unew[ix*ny+iy] =  
+                (u[(ix+1)*ny+iy] +  u[(ix-1)*ny+iy] +
+                 u[ix*ny+iy+1] +  u[ix*ny+iy-1] )/4.0;
+          
+        }
+    }
+
+//compute reduction 
+
+
+        float mydiff;
+        
+#pragma omp parallel  shared(nx,ny,u,unew, diff) private (ix,iy,mydiff)
+        {
+mydiff=0.0;
+#pragma omp for
+    for (ix = 1; ix < nx-1; ix++) {
+        for (iy = 1; iy < ny-1; iy++) {
+            if (mydiff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
+            {
+                mydiff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
+            }
+        }
+    }
+
+
+# pragma omp critical
+          {
+              if (*diff < mydiff )
+              {
+                  *diff = mydiff;
+              }
+          }
+
+
+#pragma omp for
+    for (ix = 1; ix < nx-1; ix++) {
+        for (iy = 1; iy < ny-1; iy++) {
+            u[ix*ny+iy] = unew[ix*ny+iy]; 
+        }
+    }   
+
+
+
+
+        }
+}
+
+/*****************************************************************************
+ *  Initialize Data
+ *****************************************************************************/
+void inidat(int nx, int ny, float *u, float *unew) 
+{
+    int ix, iy;
+
+    /*
+     *Set boundary data and interrior values
+     * */
+    for (ix = 0; ix < nx; ix++) 
+        for (iy = 0; iy < ny; iy++) { 
+
+            if(ix==0)
+            {
+                u[ix*ny+iy]=0.0; 
+            }
+            else          
+            if(iy==0 && ix!=0)
+            {
+                u[ix*ny+iy]=100.0;
+            }else
+           
+            if(ix==nx-1)
+            {
+                u[ix*ny+iy]=100.0;
+            }else
+
+            if(iy==ny-1 && ix!=0)
+            {   
+                u[ix*ny+iy]=100.0;
+            }else
+
+             u[ix*ny+iy]=( float ) ( 2 * nx + 2 * ny - 4 );
+        }
+}
+
+/**************************************************************************
+ * Print Data to files
+ **************************************************************************/
+void prtdat(int nx, int ny, float *u,const char *fnam)
+{
+ 
+    int ix, iy;
+    FILE *fp;
+
+    if(ITER_PRINT==0)return;
+    
+    fp = fopen(fnam, "w");
+
+    for (ix = 0 ; ix < nx; ix++) {
+        for (iy =0; iy < ny; iy++) {
+
+            fprintf(fp, "%8.3f", u[ix*ny+iy]);
+
+            if(iy!=ny-1)
+            {
+                fprintf(fp, " ");
+            }else
+            {
+                fprintf(fp, "\n");
+            }
+        }
+    }
+
+    fclose(fp);
+}
+
+
+
diff --git a/lab3/ser_heat2D.c b/lab3/ser_heat2D.c
new file mode 100644
index 0000000..bb14733
--- /dev/null
+++ b/lab3/ser_heat2D.c
@@ -0,0 +1,204 @@
+/****************************************************************************
+ * DESCRIPTION:  
+ *   Serial HEAT2D Example - C Version
+ *   This example is based on a simplified 
+ *   two-dimensional heat equation domain decomposition.  The initial 
+ *   temperature is computed to be high in the middle of the domain and 
+ *   zero at the boundaries.  The boundaries are held at zero throughout 
+ *   the simulation.  During the time-stepping, an array containing two 
+ *   domains is used; these domains alternate between old data and new data.
+ *
+ *  The physical region, and the boundary conditions, are suggested
+    by this diagram;
+
+                   u = 0
+             +------------------+
+             |                  |
+    u = 100  |                  | u = 100
+             |                  |
+             +------------------+
+                   u = 100
+
+Interrior point :
+  u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
+
+ ****************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#define N 200
+#define M 200  
+
+#define ITER_PRINT 100
+#define PRINT_DATA 1
+
+#define EPSILON 1e-1
+
+
+void update(int nx,int ny, float *u, float *unew, float * diff);
+void inidat(int nx, int ny, float *u, float *unew); 
+void prtdat(int nx, int ny, float *u,const char *fnam);
+
+
+
+
+int main(int argc, char *argv[])
+{
+
+    float diff=1.0;
+
+    float *u     = (float *)malloc(N * M * sizeof(float));
+    float *unew  = (float *)malloc(N * M * sizeof(float));
+ 
+    if(u==0 || unew ==0)
+    {
+        perror("Can't allocated data\n");
+        return -1;
+    }
+
+     printf ( "\n" );
+     printf ( "HEATED_PLATE\n" );
+     printf ( "  Serial version\n" );
+     printf ( "  A program to solve for the steady state temperature distribution\n" );
+     printf ( "  over a rectangular plate.\n" );
+     printf ( "\n" );
+     printf ( "  Spatial grid of %d by %d points.\n\n", M, N );
+
+
+    /* Initialize grid and create input file */
+    printf("Initializing grid\n");
+    
+    inidat(N, M,u,unew);
+
+    prtdat(N, M,u, "initial.dat");
+    
+   
+    printf("Start computing\n");
+
+    int iter=0;
+
+    /* 
+     *   iterate until the  new solution unew differs from the old solution u
+     *     by no more than EPSILON.
+     *     */
+     
+    while(diff> EPSILON) {
+
+        update(N, M, u, unew,&diff);
+    
+        if(iter%ITER_PRINT==0)
+        printf("Iteration %d, diff = %f\n ", iter,diff);
+
+        iter++;
+    }
+
+    prtdat(N, M, u, "final.dat");
+     
+    free(u);
+    free(unew);
+}
+
+
+
+/****************************************************************************
+ *  subroutine update
+ ****************************************************************************/
+void update(int nx,int ny, float *u, float *unew, float * diff)
+{
+    int ix, iy;
+    *diff=0.0;
+
+    for (ix = 1; ix < nx-1; ix++) {
+        for (iy = 1; iy < ny-1; iy++) {
+            unew[ix*ny+iy] =  
+                (u[(ix+1)*ny+iy] +  u[(ix-1)*ny+iy] +
+                 u[ix*ny+iy+1] +  u[ix*ny+iy-1] )/4.0;
+
+            if (*diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
+            {
+                *diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
+            }
+        }
+
+    }
+
+
+    for (ix = 1; ix < nx-1; ix++) {
+        for (iy = 1; iy < ny-1; iy++) {
+            u[ix*ny+iy] = unew[ix*ny+iy]; 
+        }
+    }   
+
+}
+
+/*****************************************************************************
+ *  Initialize Data
+ *****************************************************************************/
+void inidat(int nx, int ny, float *u, float *unew) 
+{
+    int ix, iy;
+
+    /*
+     *Set boundary data and interrior values
+     * */
+    for (ix = 0; ix < nx; ix++) 
+        for (iy = 0; iy < ny; iy++) { 
+
+            if(ix==0)
+            {
+                u[ix*ny+iy]=0.0; 
+            }
+            else          
+            if(iy==0 && ix!=0)
+            {
+                u[ix*ny+iy]=100.0;
+            }else
+           
+            if(ix==nx-1)
+            {
+                u[ix*ny+iy]=100.0;
+            }else
+
+            if(iy==ny-1 && ix!=0)
+            {   
+                u[ix*ny+iy]=100.0;
+            }else
+
+             u[ix*ny+iy]=( float ) ( 2 * nx + 2 * ny - 4 );
+        }
+}
+
+/**************************************************************************
+ * Print Data to files
+ **************************************************************************/
+void prtdat(int nx, int ny, float *u,const char *fnam)
+{
+ 
+    int ix, iy;
+    FILE *fp;
+
+    if(ITER_PRINT==0)return;
+    
+    fp = fopen(fnam, "w");
+
+    for (ix = 0 ; ix < nx; ix++) {
+        for (iy =0; iy < ny; iy++) {
+
+            fprintf(fp, "%8.3f", u[ix*ny+iy]);
+
+            if(iy!=ny-1)
+            {
+                fprintf(fp, " ");
+            }else
+            {
+                fprintf(fp, "\n");
+            }
+        }
+    }
+
+    fclose(fp);
+}
+
+
+