diff --git a/README.md b/README.md
index b60d30fe29e8b3bde5c92207ce4e1c0558d2219e..6c54f86cc5bb8f5b6be30c8fa4286ffd5a351ad9 100644
--- a/README.md
+++ b/README.md
@@ -2,3 +2,8 @@
 
 Funções e bibliotecas de uso genérico para a disciplina CI1164 (DINF/UFPR)
 
+* gnuplot: documentação e modelos de scripts para uso de gnuplot.
+* sislin: funções básicas de manipulação de matrizes e Sistemas Lineares
+* utils: funções utilitárias diversas: timestamp, geradores de
+         números pseudo-aleatórios, etc.
+
diff --git a/gnuplot/gnuplot.pdf b/gnuplot/gnuplot.pdf
index 3ff8881440946a4c547aa774e198174262999f94..8d7e815c5a47c090bc9911c7c9a41f4c7742cfdf 100644
Binary files a/gnuplot/gnuplot.pdf and b/gnuplot/gnuplot.pdf differ
diff --git a/sislin/matrix.c b/sislin/matrix.c
new file mode 100644
index 0000000000000000000000000000000000000000..e23d062506849a04525c41deb6e9168ab0c67613
--- /dev/null
+++ b/sislin/matrix.c
@@ -0,0 +1,117 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "matrix.h"
+
+/* Cria com calloc() matriz 'n x m' */
+double **criaMatriz (int n, int m)
+{
+  double **C;
+  int i, j;
+  
+  C = (double **) calloc(n,sizeof(double *));
+  C[0] = (double *) calloc(n*m,sizeof(double));
+  for (i=1; i < n; ++i)
+    C[i] = C[i-1] + m;
+
+  return C;
+}
+
+/* Retorna ponteiro para matriz inversa de matriz 'A'.
+   Matriz inversa é alocada com 'criaMatriz()'
+*/
+// Adaptação de https://www.programming-techniques.com/2011/09/numerical-methods-inverse-of-nxn-matrix.html
+double **invMatriz (double **A, int n)
+{
+  double ratio,a, **aux, **IA;
+  int i, j, k;
+
+  aux = criaMatriz(n, 2*n);
+  IA = criaMatriz(n,n);
+  
+  for(i = 0; i < n; i++) {
+    for(j = 0; j < 2*n; j++) {
+      if (j < n)
+	aux[i][j] = A[i][j];
+      else if(i==(j-n))
+	aux[i][j] = 1.0;
+      else
+	aux[i][j] = 0.0;
+    }
+  }
+
+  for(i = 0; i < n; i++){
+    for(j = 0; j < n; j++){
+      if(i!=j){
+	ratio = aux[j][i]/aux[i][i];
+	for(k = 0; k < 2*n; k++){
+	  aux[j][k] -= ratio * aux[i][k];
+	}
+      }
+    }
+  }
+
+  for(i = 0; i < n; i++){
+    a = aux[i][i];
+    for(j = n; j < 2*n; j++){
+      IA[i][j-n] = aux[i][j] / a;
+    }
+  }
+
+  freeMatriz(aux, n);
+  
+  return IA;
+
+}
+
+/* Libera matriz com 'n' linhas, alocada com 'criaMatriz()' */
+void freeMatriz (double **A, int n)
+{
+  int i;
+  
+  free(A[0]);
+  free(A);
+}
+
+/* Multiplica matrizes 'A' e 'B', ambas 'n x n' */
+double **multMatriz (double **A, double **B, int n)
+{
+  double **C;
+  int i, j,k;
+
+  C = (double **) malloc(n*sizeof(double *));
+  for (i=0; i < n; ++i) {
+    C[i] = (double *) calloc(n,sizeof(double));
+    for (j=0; j < n; ++j)
+      for (k=0; k < n; ++k)
+	C[i][j] += A[i][k] * B[k][j];
+  }
+
+  return C;
+}
+
+/* Mostra na tela matriz 'n x n' alocada com 'criaMatriz()' */
+void prnMatriz (double **A, int n)
+{
+  int i, j;
+  
+  for(i = 0; i < n; i++){
+    for(j = 0; j < n; j++){
+      printf("%15.5lg", A[i][j]);
+    }
+    printf("\n");
+  }
+}
+
+/* Preenche matriz 'n x n' alocada com 'criaMatriz()' */
+void lerMatriz (double **A, int n)
+{
+  int i, j;
+  
+  for(i = 0; i < n; i++){
+    for(j = 0; j < n; j++){
+      scanf("%lg", A[i]+j);
+    }
+  }
+}
+
diff --git a/sislin/matrix.h b/sislin/matrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..53a7bbb52dfd14b87331e5acb1617044fb88f981
--- /dev/null
+++ b/sislin/matrix.h
@@ -0,0 +1,20 @@
+/* Cria com calloc() matriz 'n x m' */
+double **criaMatriz (int n, int m);
+
+/* Retorna ponteiro para matriz inversa de matriz 'A'.
+   Matriz inversa é alocada com 'criaMatriz()'
+*/
+double **invMatriz (double **A, int n);
+
+/* Libera matriz com 'n' linhas, alocada com 'criaMatriz()' */
+void freeMatriz (double **A, int n);
+
+/* Multiplica matrizes 'A' e 'B', ambas 'n x n' */
+double **multMatriz (double **A, double **B, int n);
+
+/* Mostra na tela matriz 'n x n' alocada com 'criaMatriz()' */
+void prnMatriz (double **A, int n);
+
+/* Preenche matriz 'n x n' alocada com 'criaMatriz()' */
+void lerMatriz (double **A, int n);
+
diff --git a/sislin/sislin.c b/sislin/sislin.c
new file mode 100644
index 0000000000000000000000000000000000000000..75f649be5051939638a7a78905832b95a2d4633c
--- /dev/null
+++ b/sislin/sislin.c
@@ -0,0 +1,186 @@
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+#include "utils.h"
+#include "sislin.h"
+
+// Alocaçao de matriz em memória. 
+SistLinear_t* alocaSisLin (unsigned int n, tipoAloc_t tipo)
+{
+  SistLinear_t *SL = (SistLinear_t *) malloc(sizeof(SistLinear_t));
+  
+  if ( SL ) {
+    
+    SL->n = n;
+    SL->tipoAloc_A = tipo;
+    SL->A = (real_t **) malloc(n * sizeof(real_t *));
+    SL->b = (real_t *) malloc(n * sizeof(real_t));
+
+    if (!(SL->A) || !(SL->b)) {
+      liberaSisLin(SL);
+      return NULL;
+    }
+
+    // Matriz como vetor de N ponteiros para um único vetor com N*N elementos
+    if (tipo == pontVet) {
+      SL->A[0] = (real_t *) malloc(n * n * sizeof(real_t));
+      if (!(SL->A[0])) {
+	liberaSisLin(SL);
+	return NULL;
+      }
+	
+      for (int i=1; i < n; ++i) {
+	SL->A[i] = SL->A[i-1]+n;
+      }
+    }
+    else if (tipo == pontPont) { // Matriz  como  vetor de  N  ponteiros
+				 // para N vetores de N elementos cada
+    
+      for (int i=0; i < n; ++i)
+	SL->A[i] = (real_t *) malloc(n * sizeof(real_t));
+    }
+  }
+  
+  return (SL);
+}
+
+// Liberacao de memória
+void liberaSisLin (SistLinear_t *SL)
+{
+  if (SL) {
+    if (SL->A) {
+      if (SL->tipoAloc_A == pontVet) {
+	if (SL->A[0]) free (SL->A[0]);
+      }
+      else if (SL->tipoAloc_A == pontPont) {
+	for (int i=0; i < SL->n; ++i)
+	  free (SL->A[i]);
+      }
+      
+      free(SL->A);
+    }
+    if (SL->b) free(SL->b);
+    free(SL);
+  }
+}
+
+
+/*!
+  \brief Cria coeficientes e termos independentes do SL
+  *
+  \param SL Ponteiro para o sistema linear
+  \param tipo Tipo de sistema linear a ser criado. Pode ser: 
+     comSolucao, eqNula, eqProporcional, eqCombLinear, hilbert 
+  \param coef_max Maior valor para coeficientes e termos independentes
+*/
+void iniSisLin (SistLinear_t *SL, tipoSistLinear_t tipo, real_t coef_max)
+{
+  unsigned int n = SL->n;
+  // para gerar valores no intervalo [0,coef_max]
+  real_t invRandMax = ((real_t)coef_max / (real_t)RAND_MAX);
+
+  // inicializa vetor b
+  for (unsigned int i=0; i<n; ++i) {
+    SL->b[i] = (real_t)rand() * invRandMax;
+  }
+    
+  if (tipo == hilbert) {
+    for (unsigned int i=0; i<n; ++i) {
+      for (unsigned int j=0; j<n; ++j)  {
+	SL->A[i][j] = 1.0 / (real_t)(i+j+1);
+      }
+    }
+  }
+  else { // inicializa sistema normal e depois altera
+    // inicializa a matriz A
+    for (unsigned int i=0; i<n; ++i) {
+      for (unsigned int j=0; j<n; ++j)  {
+	SL->A[i][j] = (real_t)rand() * invRandMax;
+      }
+    }
+    if (tipo == eqNula) {
+      // sorteia eq a ser "nula"
+      unsigned int nula = rand() % n;
+      for (unsigned int j=0; j<n; ++j) {
+	SL->A[nula][j] = 0.0;
+      }
+      SL->b[nula] = 0.0;
+    } 
+    else if (tipo == eqProporcional) {
+      // sorteia eq a ser "proporcional" e valor
+      unsigned int propDst = rand() % n;
+      unsigned int propSrc = (propDst + 1) % n;
+      real_t mult = (real_t)rand() * invRandMax;
+      for (unsigned int j=0; j<n; ++j) {
+	SL->A[propDst][j] = SL->A[propSrc][j] * mult;
+      }
+      SL->b[propDst] = SL->b[propSrc] * mult;
+    } 
+    else if (tipo == eqCombLinear) {
+      // sorteia eq a ser "combLinear"
+      unsigned int combDst = rand() % n;
+      unsigned int combSrc1 = (combDst + 1) % n;
+      unsigned int combSrc2 = (combDst + 2) % n;
+      for (unsigned int j=0; j<n; ++j) {
+	SL->A[combDst][j] = SL->A[combSrc1][j] + SL->A[combSrc2][j];
+      }
+      SL->b[combDst] = SL->b[combSrc1] + SL->b[combSrc2];
+    }
+    else if (tipo == diagDominante) {
+      // aumenta o valor dos termos da diagonal principal
+      for (unsigned int i=0; i<n; ++i) {
+	real_t soma = 0.0;
+	for (unsigned int j=0; j < i; ++j) soma += SL->A[i][j];
+	for (unsigned int j=i+1; j < n; ++j) soma += SL->A[i][j];
+        SL->A[i][i] += soma;
+      }
+    }
+  }
+}
+
+
+SistLinear_t *lerSisLin (tipoAloc_t tipo)
+{
+  unsigned int n;
+  SistLinear_t *SL;
+  
+  scanf("%d",&n);
+
+  SL = alocaSisLin (n, tipo);
+  
+  for(int i=0; i < n; ++i)
+    for(int j=0; j < n; ++j)
+      scanf ("%lg", &SL->A[i][j]);
+
+  for(int i=0; i < n; ++i)
+    scanf ("%lg", &SL->b[i]);
+  
+  return SL;
+}
+
+
+void prnSisLin (SistLinear_t *SL)
+{
+  int n=SL->n;
+
+  for(int i=0; i < n; ++i) {
+    printf("\n  ");
+    for(int j=0; j < n; ++j)
+      printf ("%10g", SL->A[i][j]);
+    printf ("   |   %g", SL->b[i]);
+  }
+  printf("\n\n");
+}
+
+void prnVetor (real_t *v, unsigned int n)
+{
+  int i;
+
+  printf ("\n");
+  for(i=0; i < n; ++i)
+      printf ("%10g ", v[i]);
+  printf ("\n\n");
+
+}
+
diff --git a/sislin/sislin.h b/sislin/sislin.h
new file mode 100644
index 0000000000000000000000000000000000000000..e861beecb8040ad3b012862a9bbd5d471d727650
--- /dev/null
+++ b/sislin/sislin.h
@@ -0,0 +1,43 @@
+#ifndef __SISLIN_H__
+#define __SISLIN_H__
+
+#define COEF_MAX 32.0 // Valor máximo usado para gerar valores aleatórios de
+		      // coeficientes nos sistemas lineares.
+
+// Tipo de alocação para matrizes
+typedef enum {
+  pontPont=0, // Matriz como vetor de N ponteiros para vetores de tamanho N
+  pontVet     // Matriz como vetor de N ponteiros para um único vetor de tamanho N*N
+} tipoAloc_t;
+
+// Estrutura para definiçao de um sistema linear qualquer
+typedef struct {
+  real_t **A; // coeficientes
+  real_t *b; // termos independentes
+  unsigned int n; // tamanho do SL
+  tipoAloc_t tipoAloc_A; // tipo de alocação usada na matriz de coeficientes
+} SistLinear_t;
+
+// Tipos de matrizes de coeficientes usados pela função 'inicializaSistLinear()'
+typedef enum {
+    generico = 0,
+    hilbert,
+    diagDominante,
+    eqNula,
+    eqProporcional,
+    eqCombLinear
+} tipoSistLinear_t;
+
+
+// Alocaçao e desalocação de matrizes
+SistLinear_t* alocaSisLin (unsigned int n, tipoAloc_t tipo);
+void liberaSisLin (SistLinear_t *SL);
+void iniSisLin (SistLinear_t *SL, tipoSistLinear_t tipo, real_t coef_max);
+
+// Leitura e impressão de sistemas lineares
+SistLinear_t *lerSisLin ();
+void prnSisLin (SistLinear_t *SL);
+void prnVetor (real_t *vet, unsigned int n);
+
+#endif // __SISLIN_H__
+