---Step 1: Creating INP.File---
[ rawData,dat ] = creating_inp_file( '20x20x2.mph','mfnc.Bz',[10 13],'2mm_mfnc.Bz.inp');

function [ rawData,dat ] = creating_inp_file( COMSOL_file,Bz,sections,INP_File)

%---Loading COMSOL_file---

    model=mphload(COMSOL_file);
    dat = mpheval(model,Bz,'edim','domain','selection',sections);
    sumOfNodes=length(dat.p(1,:));
    sumOfElements=length(dat.t(1,:));
    
%---Creating inp_file---    
    
 % Number of Lines in the file
 
    iLines=36+sumOfNodes+sumOfElements+ceil(sumOfNodes/9);	

 % Cellarray for the fileinput (preallocation)
 
    rawData=cell(iLines,9);
    
 % Writting inp_file
    rawData{1,1}='*Heading';
    rawData{2,1}='**Job name: COMSOL_file';
    rawData{3,1}='**Generated by: COMSOL Version 5.1';
    rawData{4,1}='*Preprint, echo=NO, model=NO, history=NO, contact=NO';
    rawData{5,1}='*Part, name=PART-1';
    rawData{6,1}='*Node';
    
    for i=1:sumOfNodes
        rawData{6+i,1}=i;%Node numbers
        counter=1;
        for j=2:4
            rawData{6+i,j}=dat.p(counter,i);
            counter=counter+1;
        end
    end
    
    rawData{7+sumOfNodes,1}='*Element, type=C3D4, elset=set1';
    
    for i=1:sumOfElements
          rawData{7+sumOfNodes+i,1}=i;
          counter=1;
        for j=2:5
            rawData{7+sumOfNodes+i,j}=dat.t(counter,i)+1;
            counter=counter+1;
        end
    end
    
    rawData{8+sumOfNodes+sumOfElements,1}='*Nset, nest=ALL';
    counter=1;
    
    for i=1:ceil(sumOfNodes/9)
        for j=1:9
            if counter == sumOfNodes+1
                 break
            end
             rawData{8+sumOfNodes+sumOfElements+i,j}=counter;
             counter=counter+1;
        end
    end
    
     rawData{9+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*SOLID SECTION, ELSET=set1, MATERIAL=PBM';
     rawData{10+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}=1.0;
     rawData{11+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*End Part';
     rawData{12+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*Assembly, name=Assembly';
     rawData{13+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*Instance, name=PART-1-2, part=PART-1';
     rawData{14+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*End Instance';
     rawData{15+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*End Assembly';
     rawData{16+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*Material, name=PBM';
     rawData{17+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='** ----------------------------------------------------------------';
     rawData{18+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*Step, name=Step-1';
     rawData{19+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*Static';
     rawData{20+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='2., 2., 2e-05, 2.';
     rawData{21+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*Restart, write, frequency=0';
     rawData{22+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*Output, field, variable=PRESELECT';
     rawData{23+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*Output, history, variable=PRESELECT';
     rawData{24+sumOfNodes+sumOfElements+ceil(sumOfNodes/9)}='*End Step';

 % Printing inp_file
     iLines=length(rawData);
     fileID = fopen(INP_File,'w');

        for i=1:iLines
           if iscellstr(rawData(i))==1
               formatSpec='%s\n';
               fprintf(fileID,formatSpec,rawData{i});
           else
            for j=1:8
                if isempty(rawData{i,j})==0 & isempty(rawData{i,j+1})==0
                    formatSpec='%d, ';
                    fprintf(fileID,formatSpec,rawData{i,j});
                    if j==8
                        formatSpec2='%d\n';
                        fprintf(fileID,formatSpec2,rawData{i,j+1});
                    end
                else 
                    formatSpec='%d\n';
                    fprintf(fileID,formatSpec,rawData{i,j});
                    break
                end    
             end
           end
        end

      fclose(fileID);