diff --git a/code/micropython.mk b/code/micropython.mk
index d16b177..f36d1d6 100644
--- a/code/micropython.mk
+++ b/code/micropython.mk
@@ -19,6 +19,7 @@ SRC_USERMOD += $(USERMODULES_DIR)/numpy/create.c
 SRC_USERMOD += $(USERMODULES_DIR)/numpy/fft/fft.c
 SRC_USERMOD += $(USERMODULES_DIR)/numpy/fft/fft_tools.c
 SRC_USERMOD += $(USERMODULES_DIR)/numpy/filter.c
+SRC_USERMOD += $(USERMODULES_DIR)/numpy/io/io.c
 SRC_USERMOD += $(USERMODULES_DIR)/numpy/linalg/linalg.c
 SRC_USERMOD += $(USERMODULES_DIR)/numpy/linalg/linalg_tools.c
 SRC_USERMOD += $(USERMODULES_DIR)/numpy/numerical.c
diff --git a/code/numpy/io/io.c b/code/numpy/io/io.c
new file mode 100644
index 0000000..4beefa1
--- /dev/null
+++ b/code/numpy/io/io.c
@@ -0,0 +1,375 @@
+/*
+ * This file is part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2022 Zoltán Vörös
+*/
+
+#include <string.h>
+#include <stdio.h>
+
+#include "py/builtin.h"
+#include "py/obj.h"
+#include "py/runtime.h"
+#include "py/stream.h"
+
+#include "../../ndarray.h"
+#include "io.h"
+
+#define ULAB_IO_BUFFER_SIZE         128
+
+#define ULAB_IO_NULL_ENDIAN         0
+#define ULAB_IO_LITTLE_ENDIAN       1
+#define ULAB_IO_BIG_ENDIAN          2
+
+#if ULAB_NUMPY_HAS_LOAD
+static void io_read_(mp_obj_t stream, const mp_stream_p_t *stream_p, char *buffer, char *string, uint16_t len, int *error) {
+    size_t read = stream_p->read(stream, buffer, len, error);
+    bool fail = false;
+    if(read == len) {
+        if(string != NULL) {
+            if(memcmp(buffer, string, len) != 0) {
+                fail = true;
+            }
+        }
+    } else {
+        fail = true;
+    }
+    if(fail) {
+        stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, error);
+        mp_raise_msg(&mp_type_RuntimeError, translate("corrupted file"));
+    }
+}
+
+static mp_obj_t io_load(mp_obj_t file) {
+    if(!mp_obj_is_str(file)) {
+        mp_raise_TypeError(translate("wrong input type"));
+    }
+
+    int error;
+    char *buffer = m_new(char, ULAB_IO_BUFFER_SIZE);
+
+    // test for endianness
+    uint16_t x = 1;
+    int8_t native_endianness = (x >> 8) == 1 ? ULAB_IO_BIG_ENDIAN : ULAB_IO_LITTLE_ENDIAN;
+
+    mp_obj_t open_args[2] = {
+        file,
+        MP_OBJ_NEW_QSTR(MP_QSTR_rb)
+    };
+
+    mp_obj_t stream = mp_builtin_open(2, open_args, (mp_map_t *)&mp_const_empty_map);
+    const mp_stream_p_t *stream_p = mp_get_stream(stream);
+
+    // read header
+    // magic string
+    io_read_(stream, stream_p, buffer, "\x93NUMPY", 6, &error);
+    // simply discard the version number
+    io_read_(stream, stream_p, buffer, NULL, 2, &error);
+    // header length, represented as a little endian uint16 (0x76, 0x00)
+    io_read_(stream, stream_p, buffer, NULL, 2, &error);
+
+    uint16_t header_length = buffer[1];
+    header_length <<= 8;
+    header_length += buffer[0];
+
+    // beginning of the dictionary describing the array
+    io_read_(stream, stream_p, buffer, "{'descr': '", 11, &error);
+    uint8_t dtype;
+
+    io_read_(stream, stream_p, buffer, NULL, 1, &error);
+    uint8_t endianness = ULAB_IO_NULL_ENDIAN;
+    if(*buffer == '<') {
+        endianness = ULAB_IO_LITTLE_ENDIAN;
+    } else if(*buffer == '>') {
+        endianness = ULAB_IO_BIG_ENDIAN;
+    }
+
+    io_read_(stream, stream_p, buffer, NULL, 2, &error);
+    if(memcmp(buffer, "u1", 2) == 0) {
+        dtype = NDARRAY_UINT8;
+    } else if(memcmp(buffer, "i1", 2) == 0) {
+        dtype = NDARRAY_INT8;
+    } else if(memcmp(buffer, "u2", 2) == 0) {
+        dtype = NDARRAY_UINT16;
+    } else if(memcmp(buffer, "i2", 2) == 0) {
+        dtype = NDARRAY_INT16;
+    }
+    #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
+    else if(memcmp(buffer, "f4", 2) == 0) {
+        dtype = NDARRAY_FLOAT;
+    }
+    #else
+    else if(memcmp(buffer, "f8", 2) == 0) {
+        dtype = NDARRAY_FLOAT;
+    }
+    #endif
+    #if ULAB_SUPPORTS_COMPLEX
+    #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
+    else if(memcmp(buffer, "c8", 2) == 0) {
+        dtype = NDARRAY_COMPLEX;
+    }
+    #else
+    else if(memcmp(buffer, "c16", 3) == 0) {
+        dtype = NDARRAY_COMPLEX;
+    }
+    #endif
+    #endif /* ULAB_SUPPORT_COPMLEX */
+    else {
+        stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
+        mp_raise_TypeError(translate("wrong dtype"));
+    }
+
+    io_read_(stream, stream_p, buffer, "', 'fortran_order': False, 'shape': (", 37, &error);
+
+    size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
+    memset(shape, 0, sizeof(size_t) * ULAB_MAX_DIMS);
+
+    uint16_t bytes_to_read = MIN(ULAB_IO_BUFFER_SIZE, header_length - 51);
+    // bytes_to_read is 128 at most. This should be enough to contain a
+    // maximum of 4 size_t numbers plus the delimiters
+    io_read_(stream, stream_p, buffer, NULL, bytes_to_read, &error);
+    char *needle = buffer;
+    uint8_t ndim = 0;
+
+    // find out the number of dimensions by counting the commas in the string
+    while(1) {
+        if(*needle == ',') {
+            ndim++;
+            if(needle[1] == ')') {
+                break;
+            }
+        } else if((*needle == ')') && (ndim > 0)) {
+            ndim++;
+            break;
+        }
+        needle++;
+    }
+
+    needle = buffer;
+    for(uint8_t i = 0; i < ndim; i++) {
+        size_t number = 0;
+        // trivial number parsing here
+        while(1) {
+            if((*needle == ' ') || (*needle == '\t')) {
+                needle++;
+            }
+            if((*needle > 47) && (*needle < 58)) {
+                number = number * 10 + (*needle - 48);
+            } else if((*needle == ',') || (*needle == ')')) {
+                break;
+            }
+            else {
+                stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
+                mp_raise_msg(&mp_type_RuntimeError, translate("corrupted file"));
+            }
+            needle++;
+        }
+        needle++;
+        shape[ULAB_MAX_DIMS - ndim + i] = number;
+    }
+
+    // strip the rest of the header
+    if((bytes_to_read + 51) < header_length) {
+        io_read_(stream, stream_p, buffer, NULL, header_length - (bytes_to_read + 51), &error);
+    }
+
+    ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(ndim, shape, dtype);
+    char *array = (char *)ndarray->array;
+
+    size_t read = stream_p->read(stream, array, ndarray->len * ndarray->itemsize, &error);
+    if(read != ndarray->len * ndarray->itemsize) {
+        stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
+        mp_raise_msg(&mp_type_RuntimeError, translate("corrupted file"));
+    }
+
+    stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
+    m_del(char, buffer, ULAB_IO_BUFFER_SIZE);
+
+    // swap the bytes, if necessary
+    if((native_endianness != endianness) && (dtype != NDARRAY_UINT8) && (dtype != NDARRAY_INT8)) {
+        uint8_t sz = ndarray->itemsize;
+        char *tmpbuff = NULL;
+
+        #if ULAB_SUPPORTS_COMPLEX
+        if(dtype == NDARRAY_COMPLEX) {
+            // work with the floating point real and imaginary parts
+            sz /= 2;
+            tmpbuff = m_new(char, sz);
+            for(size_t i = 0; i < ndarray->len; i++) {
+                for(uint8_t k = 0; k < 2; k++) {
+                    tmpbuff += sz;
+                    for(uint8_t j = 0; j < sz; j++) {
+                        memcpy(--tmpbuff, array++, 1);
+                    }
+                    memcpy(array-sz, tmpbuff, sz);
+                }
+            }
+        } else {
+        #endif
+            tmpbuff = m_new(char, sz);
+            for(size_t i = 0; i < ndarray->len; i++) {
+                tmpbuff += sz;
+                for(uint8_t j = 0; j < sz; j++) {
+                    memcpy(--tmpbuff, array++, 1);
+                }
+                memcpy(array-sz, tmpbuff, sz);
+            }
+        #if ULAB_SUPPORTS_COMPLEX
+        }
+        #endif
+        m_del(char, tmpbuff, sz);
+    }
+
+    return MP_OBJ_FROM_PTR(ndarray);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_1(io_load_obj, io_load);
+#endif /* ULAB_NUMPY_HAS_LOAD */
+
+#if ULAB_NUMPY_HAS_SAVE
+static mp_obj_t io_save(mp_obj_t file, mp_obj_t ndarray_) {
+    if(!mp_obj_is_str(file) || !mp_obj_is_type(ndarray_, &ulab_ndarray_type)) {
+        mp_raise_TypeError(translate("wrong input type"));
+    }
+
+    ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(ndarray_);
+    int error;
+    char *buffer = m_new(char, ULAB_IO_BUFFER_SIZE);
+    uint8_t offset = 0;
+
+    // test for endianness
+    uint16_t x = 1;
+    int8_t native_endiannes = (x >> 8) == 1 ? '>' : '<';
+
+    mp_obj_t open_args[2] = {
+        file,
+        MP_OBJ_NEW_QSTR(MP_QSTR_wb)
+    };
+
+    mp_obj_t stream = mp_builtin_open(2, open_args, (mp_map_t *)&mp_const_empty_map);
+    const mp_stream_p_t *stream_p = mp_get_stream(stream);
+
+    // write header;
+    // magic string + header length, which is always 128 - 10 = 118, represented as a little endian uint16 (0x76, 0x00)
+    // + beginning of the dictionary describing the array
+    memcpy(buffer, "\x93NUMPY\x01\x00\x76\x00{'descr': '", 21);
+    offset += 21;
+
+    buffer[offset] = native_endiannes;
+    if((ndarray->dtype == NDARRAY_UINT8) || (ndarray->dtype == NDARRAY_INT8)) {
+        // for single-byte data, the endianness doesn't matter
+        buffer[offset] = '|';
+    }
+    offset++;
+    switch(ndarray->dtype) {
+        case NDARRAY_UINT8:
+            memcpy(buffer+offset, "u1", 2);
+            break;
+        case NDARRAY_INT8:
+            memcpy(buffer+offset, "i1", 2);
+            break;
+        case NDARRAY_UINT16:
+            memcpy(buffer+offset, "u2", 2);
+            break;
+        case NDARRAY_INT16:
+            memcpy(buffer+offset, "i2", 2);
+            break;
+        case NDARRAY_FLOAT:
+            #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
+            memcpy(buffer+offset, "f4", 2);
+            #else
+            memcpy(buffer+offset, "f8", 2);
+            #endif
+            break;
+        #if ULAB_SUPPORTS_COMPLEX
+        case NDARRAY_COMPLEX:
+            #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
+            memcpy(buffer+offset, "c8", 2);
+            #else
+            memcpy(buffer+offset, "c16", 3);
+            offset++;
+            #endif
+            break;
+        #endif
+    }
+
+    offset += 2;
+    memcpy(buffer+offset, "', 'fortran_order': False, 'shape': (", 37);
+    offset += 37;
+
+    if(ndarray->ndim == 1) {
+        offset += sprintf(buffer+offset, "%zu,", ndarray->shape[ULAB_MAX_DIMS - 1]);
+    } else {
+        for(uint8_t i = ndarray->ndim; i > 1; i--) {
+            offset += sprintf(buffer+offset, "%zu, ", ndarray->shape[ULAB_MAX_DIMS - i]);
+        }
+        offset += sprintf(buffer+offset, "%zu", ndarray->shape[ULAB_MAX_DIMS - 1]);
+    }
+    memcpy(buffer+offset, "), }", 4);
+    offset += 4;
+    // pad with space till the very end
+    memset(buffer+offset, 32, ULAB_IO_BUFFER_SIZE - offset - 1);
+    buffer[ULAB_IO_BUFFER_SIZE - 1] = '\n';
+    stream_p->write(stream, buffer, ULAB_IO_BUFFER_SIZE, &error);
+
+    // write the array data
+    uint8_t sz = ndarray->itemsize;
+    offset = 0;
+
+    uint8_t *array = (uint8_t *)ndarray->array;
+
+    #if ULAB_MAX_DIMS > 3
+    size_t i = 0;
+    do {
+    #endif
+        #if ULAB_MAX_DIMS > 2
+        size_t j = 0;
+        do {
+        #endif
+            #if ULAB_MAX_DIMS > 1
+            size_t k = 0;
+            do {
+            #endif
+                size_t l = 0;
+                do {
+                    memcpy(buffer+offset, array, sz);
+                    offset += sz;
+                    if(offset == ULAB_IO_BUFFER_SIZE) {
+                        stream_p->write(stream, buffer, offset, &error);
+                        offset = 0;
+                    }
+                    array += ndarray->strides[ULAB_MAX_DIMS - 1];
+                    l++;
+                } while(l <  ndarray->shape[ULAB_MAX_DIMS - 1]);
+            #if ULAB_MAX_DIMS > 1
+                array -= ndarray->strides[ULAB_MAX_DIMS - 1] * ndarray->shape[ULAB_MAX_DIMS-1];
+                array += ndarray->strides[ULAB_MAX_DIMS - 2];
+                k++;
+            } while(k <  ndarray->shape[ULAB_MAX_DIMS - 2]);
+            #endif
+        #if ULAB_MAX_DIMS > 2
+            array -= ndarray->strides[ULAB_MAX_DIMS - 2] * ndarray->shape[ULAB_MAX_DIMS-2];
+            array += ndarray->strides[ULAB_MAX_DIMS - 3];
+            j++;
+        } while(j <  ndarray->shape[ULAB_MAX_DIMS - 3]);
+        #endif
+    #if ULAB_MAX_DIMS > 3
+        array -= ndarray->strides[ULAB_MAX_DIMS - 3] * ndarray->shape[ULAB_MAX_DIMS-3];
+        array += ndarray->strides[ULAB_MAX_DIMS - 4];
+        i++;
+    } while(i <  ndarray->shape[ULAB_MAX_DIMS - 4]);
+    #endif
+
+    stream_p->write(stream, buffer, offset, &error);
+    stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
+
+    m_del(char, buffer, ULAB_IO_BUFFER_SIZE);
+    return mp_const_none;
+}
+
+MP_DEFINE_CONST_FUN_OBJ_2(io_save_obj, io_save);
+#endif /* ULAB_NUMPY_HAS_SAVE */
diff --git a/code/numpy/io/io.h b/code/numpy/io/io.h
new file mode 100644
index 0000000..a9dcdfc
--- /dev/null
+++ b/code/numpy/io/io.h
@@ -0,0 +1,17 @@
+/*
+ * This file is part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2022 Zoltán Vörös
+*/
+
+#ifndef _ULAB_IO_
+#define _ULAB_IO_
+
+MP_DECLARE_CONST_FUN_OBJ_2(io_save_obj);
+MP_DECLARE_CONST_FUN_OBJ_1(io_load_obj);
+
+#endif
\ No newline at end of file
diff --git a/code/numpy/numpy.c b/code/numpy/numpy.c
index 25600db..1e593b5 100644
--- a/code/numpy/numpy.c
+++ b/code/numpy/numpy.c
@@ -23,6 +23,7 @@
 #include "create.h"
 #include "fft/fft.h"
 #include "filter.h"
+#include "io/io.h"
 #include "linalg/linalg.h"
 #include "numerical.h"
 #include "stats.h"
@@ -256,6 +257,9 @@ static const mp_rom_map_elem_t ulab_numpy_globals_table[] = {
     #if ULAB_NUMPY_HAS_FLIP
         { MP_OBJ_NEW_QSTR(MP_QSTR_flip), (mp_obj_t)&numerical_flip_obj },
     #endif
+    #if ULAB_NUMPY_HAS_LOAD
+        { MP_OBJ_NEW_QSTR(MP_QSTR_load), (mp_obj_t)&io_load_obj },
+    #endif
     #if ULAB_NUMPY_HAS_MINMAX
         { MP_OBJ_NEW_QSTR(MP_QSTR_max), (mp_obj_t)&numerical_max_obj },
     #endif
@@ -271,6 +275,9 @@ static const mp_rom_map_elem_t ulab_numpy_globals_table[] = {
     #if ULAB_NUMPY_HAS_ROLL
         { MP_OBJ_NEW_QSTR(MP_QSTR_roll), (mp_obj_t)&numerical_roll_obj },
     #endif
+    #if ULAB_NUMPY_HAS_SAVE
+        { MP_OBJ_NEW_QSTR(MP_QSTR_save), (mp_obj_t)&io_save_obj },
+    #endif
     #if ULAB_NUMPY_HAS_SIZE
         { MP_OBJ_NEW_QSTR(MP_QSTR_size), (mp_obj_t)&transform_size_obj },
     #endif
diff --git a/code/ulab.c b/code/ulab.c
index 4d0f490..8d4ae70 100644
--- a/code/ulab.c
+++ b/code/ulab.c
@@ -33,7 +33,7 @@
 #include "user/user.h"
 #include "utils/utils.h"
 
-#define ULAB_VERSION 4.2.1
+#define ULAB_VERSION 4.3.0
 #define xstr(s) str(s)
 #define str(s) #s
 
diff --git a/code/ulab.h b/code/ulab.h
index d5750bc..a782da2 100644
--- a/code/ulab.h
+++ b/code/ulab.h
@@ -462,6 +462,10 @@
 #define ULAB_NUMPY_HAS_INTERP           (1)
 #endif
 
+#ifndef ULAB_NUMPY_HAS_LOAD
+#define ULAB_NUMPY_HAS_LOAD             (1)
+#endif
+
 #ifndef ULAB_NUMPY_HAS_MEAN
 #define ULAB_NUMPY_HAS_MEAN             (1)
 #endif
@@ -486,6 +490,10 @@
 #define ULAB_NUMPY_HAS_ROLL             (1)
 #endif
 
+#ifndef ULAB_NUMPY_HAS_SAVE
+#define ULAB_NUMPY_HAS_SAVE             (1)
+#endif
+
 #ifndef ULAB_NUMPY_HAS_SIZE
 #define ULAB_NUMPY_HAS_SIZE             (1)
 #endif
diff --git a/docs/manual/source/conf.py b/docs/manual/source/conf.py
index 1275760..44f0134 100644
--- a/docs/manual/source/conf.py
+++ b/docs/manual/source/conf.py
@@ -27,7 +27,7 @@ copyright = '2019-2022, Zoltán Vörös and contributors'
 author = 'Zoltán Vörös'
 
 # The full version, including alpha/beta/rc tags
-release = '4.2.0'
+release = '4.3.0'
 
 
 # -- General configuration ---------------------------------------------------
diff --git a/docs/manual/source/numpy-functions.rst b/docs/manual/source/numpy-functions.rst
index a18e12c..a9a8111 100644
--- a/docs/manual/source/numpy-functions.rst
+++ b/docs/manual/source/numpy-functions.rst
@@ -11,7 +11,7 @@ the firmware was compiled with complex support.
 3.  `numpy.argmax <#argmax>`__
 4.  `numpy.argmin <#argmin>`__
 5.  `numpy.argsort <#argsort>`__
-6.  `numpy.asarray <#asarray>`__
+6.  `numpy.asarray\* <#asarray>`__
 7.  `numpy.clip <#clip>`__
 8.  `numpy.compress\* <#compress>`__
 9.  `numpy.conjugate\* <#conjugate>`__
@@ -25,25 +25,27 @@ the firmware was compiled with complex support.
 17. `numpy.interp <#interp>`__
 18. `numpy.isfinite <#isfinite>`__
 19. `numpy.isinf <#isinf>`__
-20. `numpy.max <#max>`__
-21. `numpy.maximum <#maximum>`__
-22. `numpy.mean <#mean>`__
-23. `numpy.median <#median>`__
-24. `numpy.min <#min>`__
-25. `numpy.minimum <#minimum>`__
-26. `numpy.not_equal <#equal>`__
-27. `numpy.polyfit <#polyfit>`__
-28. `numpy.polyval <#polyval>`__
-29. `numpy.real\* <#real>`__
-30. `numpy.roll <#roll>`__
-31. `numpy.size <#size>`__
-32. `numpy.sort <#sort>`__
-33. `numpy.sort_complex\* <#sort_complex>`__
-34. `numpy.std <#std>`__
-35. `numpy.sum <#sum>`__
-36. `numpy.trace <#trace>`__
-37. `numpy.trapz <#trapz>`__
-38. `numpy.where <#where>`__
+20. `numpy.load <#load>`__
+21. `numpy.max <#max>`__
+22. `numpy.maximum <#maximum>`__
+23. `numpy.mean <#mean>`__
+24. `numpy.median <#median>`__
+25. `numpy.min <#min>`__
+26. `numpy.minimum <#minimum>`__
+27. `numpy.not_equal <#equal>`__
+28. `numpy.polyfit <#polyfit>`__
+29. `numpy.polyval <#polyval>`__
+30. `numpy.real\* <#real>`__
+31. `numpy.roll <#roll>`__
+32. `numpy.save <#save>`__
+33. `numpy.size <#size>`__
+34. `numpy.sort <#sort>`__
+35. `numpy.sort_complex\* <#sort_complex>`__
+36. `numpy.std <#std>`__
+37. `numpy.sum <#sum>`__
+38. `numpy.trace <#trace>`__
+39. `numpy.trapz <#trapz>`__
+40. `numpy.where <#where>`__
 
 all
 ---
@@ -982,6 +984,39 @@ positions, where the input is infinite. Integer types return the
     
 
 
+load
+----
+
+``numpy``:
+https://docs.scipy.org/doc/numpy/reference/generated/numpy.load.html
+
+The function reads data from a file in ``numpy``\ ’s
+`platform-independent
+format <https://numpy.org/doc/stable/reference/generated/numpy.lib.format.html#module-numpy.lib.format>`__,
+and returns the generated array. If the endianness of the data in the
+file and the microcontroller differ, the bytes are automatically
+swapped.
+
+.. code::
+        
+    # code to be run in micropython
+    
+    from ulab import numpy as np
+    
+    a = np.load('a.npy')
+    print(a)
+
+.. parsed-literal::
+
+    array([[0.0, 1.0, 2.0, 3.0, 4.0],
+           [5.0, 6.0, 7.0, 8.0, 9.0],
+           [10.0, 11.0, 12.0, 13.0, 14.0],
+           [15.0, 16.0, 17.0, 18.0, 19.0],
+           [20.0, 21.0, 22.0, 23.0, 24.0]], dtype=float64)
+    
+    
+
+
 mean
 ----
 
@@ -1430,6 +1465,25 @@ Vertical rolls require two internal copies of single columns.
     
 
 
+save
+----
+
+``numpy``:
+https://docs.scipy.org/doc/numpy/reference/generated/numpy.save.html
+
+With the help of this function, numerical array can be save in
+``numpy``\ ’s `platform-independent
+format <https://numpy.org/doc/stable/reference/generated/numpy.lib.format.html#module-numpy.lib.format>`__.
+
+The function takes two positional arguments, the name of the output
+file, and the array.
+
+.. code::
+
+    # code to be run in CPython
+    
+    a = np.array(range(25)).reshape((5, 5))
+    np.save('a.npy', a)
 size
 ----
 
diff --git a/docs/numpy-functions.ipynb b/docs/numpy-functions.ipynb
index 815e4be..92660e2 100644
--- a/docs/numpy-functions.ipynb
+++ b/docs/numpy-functions.ipynb
@@ -34,8 +34,8 @@
    "execution_count": 1,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-01-14T19:55:15.200755Z",
-     "start_time": "2022-01-14T19:55:15.193656Z"
+     "end_time": "2022-01-15T08:50:03.152522Z",
+     "start_time": "2022-01-15T08:50:03.141317Z"
     }
    },
    "outputs": [],
@@ -52,8 +52,8 @@
    "execution_count": 2,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-01-14T19:55:17.871864Z",
-     "start_time": "2022-01-14T19:55:17.858935Z"
+     "end_time": "2022-01-15T08:50:04.183008Z",
+     "start_time": "2022-01-15T08:50:04.162758Z"
     }
    },
    "outputs": [],
@@ -237,7 +237,7 @@
     "1. [numpy.argmax](#argmax)\n",
     "1. [numpy.argmin](#argmin)\n",
     "1. [numpy.argsort](#argsort)\n",
-    "1. [numpy.asarray](#asarray)\n",
+    "1. [numpy.asarray*](#asarray)\n",
     "1. [numpy.clip](#clip)\n",
     "1. [numpy.compress*](#compress)\n",
     "1. [numpy.conjugate*](#conjugate)\n",
@@ -251,6 +251,7 @@
     "1. [numpy.interp](#interp)\n",
     "1. [numpy.isfinite](#isfinite)\n",
     "1. [numpy.isinf](#isinf)\n",
+    "1. [numpy.load](#load)\n",
     "1. [numpy.max](#max)\n",
     "1. [numpy.maximum](#maximum)\n",
     "1. [numpy.mean](#mean)\n",
@@ -262,6 +263,7 @@
     "1. [numpy.polyval](#polyval)\n",
     "1. [numpy.real*](#real)\n",
     "1. [numpy.roll](#roll)\n",
+    "1. [numpy.save](#save)\n",
     "1. [numpy.size](#size)\n",
     "1. [numpy.sort](#sort)\n",
     "1. [numpy.sort_complex*](#sort_complex)\n",
@@ -1418,6 +1420,50 @@
     "print('\\nisinf(c):\\n', np.isinf(c))"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## load\n",
+    "\n",
+    "`numpy`: https://docs.scipy.org/doc/numpy/reference/generated/numpy.load.html\n",
+    "\n",
+    "The function reads data from a file in `numpy`'s [platform-independent format](https://numpy.org/doc/stable/reference/generated/numpy.lib.format.html#module-numpy.lib.format), and returns the generated array. If the endianness of the data in the file and the microcontroller differ, the bytes are automatically swapped."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-01-12T19:11:10.361592Z",
+     "start_time": "2022-01-12T19:11:10.342439Z"
+    }
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "array([[0.0, 1.0, 2.0, 3.0, 4.0],\n",
+      "       [5.0, 6.0, 7.0, 8.0, 9.0],\n",
+      "       [10.0, 11.0, 12.0, 13.0, 14.0],\n",
+      "       [15.0, 16.0, 17.0, 18.0, 19.0],\n",
+      "       [20.0, 21.0, 22.0, 23.0, 24.0]], dtype=float64)\n",
+      "\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "%%micropython -unix 1\n",
+    "\n",
+    "from ulab import numpy as np\n",
+    "\n",
+    "a = np.load('a.npy')\n",
+    "print(a)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -1999,6 +2045,34 @@
     "print(\"\\na rolled with None:\\n\", a)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## save\n",
+    "\n",
+    "`numpy`: https://docs.scipy.org/doc/numpy/reference/generated/numpy.save.html\n",
+    "\n",
+    "With the help of this function, numerical array can be save in `numpy`'s [platform-independent format](https://numpy.org/doc/stable/reference/generated/numpy.lib.format.html#module-numpy.lib.format).\n",
+    "\n",
+    "The function takes two positional arguments, the name of the output file, and the array. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-01-15T08:51:08.827144Z",
+     "start_time": "2022-01-15T08:51:08.813813Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "a = np.array(range(25)).reshape((5, 5))\n",
+    "np.save('a.npy', a)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -2013,8 +2087,8 @@
    "execution_count": 3,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-01-14T19:58:44.044501Z",
-     "start_time": "2022-01-14T19:58:44.034585Z"
+     "end_time": "2022-01-15T08:50:57.254168Z",
+     "start_time": "2022-01-15T08:50:57.245772Z"
     }
    },
    "outputs": [
diff --git a/docs/ulab-change-log.md b/docs/ulab-change-log.md
index 9f13aab..e85b407 100644
--- a/docs/ulab-change-log.md
+++ b/docs/ulab-change-log.md
@@ -1,3 +1,9 @@
+Wed, 19 Jan 2022
+
+    version 4.3.0
+
+    implement numpy.save, numpy.load
+
 Tue, 18 Jan 2022
 
 version 4.2.1
@@ -10,6 +16,12 @@ version 4.2.0
 
     add numpy.size, asarray
 
+Wed, 12 Jan 2022
+
+    version 4.2.0
+
+    implement numpy.save, numpy.load
+
 Wed, 12 Jan 2022
 
 version 4.1.1
diff --git a/docs/ulab-convert.ipynb b/docs/ulab-convert.ipynb
index 3da26ed..389ec55 100644
--- a/docs/ulab-convert.ipynb
+++ b/docs/ulab-convert.ipynb
@@ -17,8 +17,8 @@
    "execution_count": 1,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-01-14T19:54:52.171096Z",
-     "start_time": "2022-01-14T19:54:52.162815Z"
+     "end_time": "2022-01-15T08:48:23.883953Z",
+     "start_time": "2022-01-15T08:48:23.877040Z"
     }
    },
    "outputs": [
@@ -61,7 +61,7 @@
     "author = 'Zoltán Vörös'\n",
     "\n",
     "# The full version, including alpha/beta/rc tags\n",
-    "release = '4.2.0'\n",
+    "release = '4.3.0'\n",
     "\n",
     "\n",
     "# -- General configuration ---------------------------------------------------\n",
@@ -218,8 +218,8 @@
    "execution_count": 2,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-01-14T20:05:37.425494Z",
-     "start_time": "2022-01-14T20:05:35.620545Z"
+     "end_time": "2022-01-15T08:48:32.207113Z",
+     "start_time": "2022-01-15T08:48:32.051714Z"
     }
    },
    "outputs": [],
@@ -256,11 +256,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 3,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-01-14T20:06:04.832792Z",
-     "start_time": "2022-01-14T20:06:00.259738Z"
+     "end_time": "2022-01-15T08:52:20.686225Z",
+     "start_time": "2022-01-15T08:52:16.125014Z"
     }
    },
    "outputs": [],
diff --git a/tests/2d/numpy/load_save.py b/tests/2d/numpy/load_save.py
new file mode 100644
index 0000000..6fb9d2a
--- /dev/null
+++ b/tests/2d/numpy/load_save.py
@@ -0,0 +1,14 @@
+try:
+    from ulab import numpy as np
+except:
+    import numpy as np
+
+dtypes = (np.uint8, np.int8, np.uint16, np.int16, np.float)
+
+for dtype in dtypes:
+    a = np.array(range(25), dtype=dtype)
+    b = a.reshape((5, 5))
+    np.save('out.npy', a)
+    print(np.load('out.npy'))
+    np.save('out.npy', b)
+    print(np.load('out.npy'))
diff --git a/tests/2d/numpy/load_save.py.exp b/tests/2d/numpy/load_save.py.exp
new file mode 100644
index 0000000..71ca601
--- /dev/null
+++ b/tests/2d/numpy/load_save.py.exp
@@ -0,0 +1,30 @@
+array([0, 1, 2, ..., 22, 23, 24], dtype=uint8)
+array([[0, 1, 2, 3, 4],
+       [5, 6, 7, 8, 9],
+       [10, 11, 12, 13, 14],
+       [15, 16, 17, 18, 19],
+       [20, 21, 22, 23, 24]], dtype=uint8)
+array([0, 1, 2, ..., 22, 23, 24], dtype=int8)
+array([[0, 1, 2, 3, 4],
+       [5, 6, 7, 8, 9],
+       [10, 11, 12, 13, 14],
+       [15, 16, 17, 18, 19],
+       [20, 21, 22, 23, 24]], dtype=int8)
+array([0, 1, 2, ..., 22, 23, 24], dtype=uint16)
+array([[0, 1, 2, 3, 4],
+       [5, 6, 7, 8, 9],
+       [10, 11, 12, 13, 14],
+       [15, 16, 17, 18, 19],
+       [20, 21, 22, 23, 24]], dtype=uint16)
+array([0, 1, 2, ..., 22, 23, 24], dtype=int16)
+array([[0, 1, 2, 3, 4],
+       [5, 6, 7, 8, 9],
+       [10, 11, 12, 13, 14],
+       [15, 16, 17, 18, 19],
+       [20, 21, 22, 23, 24]], dtype=int16)
+array([0.0, 1.0, 2.0, ..., 22.0, 23.0, 24.0], dtype=float64)
+array([[0.0, 1.0, 2.0, 3.0, 4.0],
+       [5.0, 6.0, 7.0, 8.0, 9.0],
+       [10.0, 11.0, 12.0, 13.0, 14.0],
+       [15.0, 16.0, 17.0, 18.0, 19.0],
+       [20.0, 21.0, 22.0, 23.0, 24.0]], dtype=float64)